{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Jaakkola & Millner, Non-dogmatic climate policy, JAERE.\n",
    "# Code runs on Julia 1.0.5\n",
    "\n",
    "# The code is not fully automatic and requires the analyst to manually run it to generate the desired output.\n",
    "# The code is separated out into various 'tasks' (the beginning of each new task is indicated in the comments)\n",
    "#   and each task contains multiple blocks.  Intermediate saving or loading of the data is sometimes useful.\n",
    "#   Instructions on how to solve each task, which blocks to skip etc. are indicated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Packages: All packages used are built using the blocks at the very end of the code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Task: initialise the model\n",
    "\n",
    "# Load the model run parameters based on the period length and the number of agent types\n",
    "# This initialises the model for the main parameters and non-dogmatic runs.\n",
    "\n",
    "using DataFrames, CSV;\n",
    "\n",
    "# The model run type is set up using the following parameters.\n",
    "#   A subfolder named according to the convention 'Y10T250N10' should be set up in advance\n",
    "#     and contain the parameter file 'preferences.csv', which contains the type preference parameters and the switching probability matrix.\n",
    "#     We provide the parameters and probabilities used in the main calibration, as well as the Matlab file used to generate them.\n",
    "#   Y stands for period length, T for time horizon, N for number of agents\n",
    "#   Additional subfolders with suffix \"robust\" required with the preferences computed using Mahalanobis weights\n",
    "period_length = 10;\n",
    "agent_type_n = 10;\n",
    "time_horizon = 500;\n",
    "\n",
    "# Set to \"robust\" for Mahalanobis weights.\n",
    "#model_run_string = \"robust\"\n",
    "model_run_string = \"\"\n",
    "\n",
    "type_labels = Vector{String}(undef, agent_type_n + 2)\n",
    "for i in eachindex(type_labels)\n",
    "    type_labels[i] = string(\"f\", i-2)\n",
    "end\n",
    "type_labels[1] = \"beta\"\n",
    "type_labels[2] = \"eta\"\n",
    "\n",
    "model_instance = string(\"Y$(time_horizon)T$(period_length)N$(agent_type_n)\", model_run_string);\n",
    "typetable = CSV.read(string(\".\\\\\", model_instance, \"\\\\preferences.csv\"), DataFrame, header = type_labels)\n",
    "\n",
    "dogmatism_runs = div(size(typetable)[1],(agent_type_n + 1))\n",
    "\n",
    "modelrun = NamedTuple[(\n",
    "                degree = typetable[(agent_type_n + 1)*(i-1)+1,1],\n",
    "                βs = typetable[(agent_type_n + 1)*(i-1)+2:(agent_type_n + 1)*(i-1)+agent_type_n+1, 1],\n",
    "                ηs = typetable[(agent_type_n + 1)*(i-1)+2:(agent_type_n + 1)*(i-1)+agent_type_n+1, 2],\n",
    "                weights = convert(Matrix, typetable[(agent_type_n + 1)*(i-1)+2:(agent_type_n + 1)*(i-1)+agent_type_n+1, 3:agent_type_n+2]))\n",
    "                for i in 1:dogmatism_runs]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Alternative initialisation for Stern and Nordhaus\n",
    "#   This is used to initialise a model which will run dogmatic 'Nordhaus' and 'Stern'\n",
    "#   scenarios (types 1 and 2, respectively)\n",
    "# RUN THIS INSTEAD OF THE PREVIOUS BLOCK IF GENERATING THE NORDHAUS/STERN PATHS, OTHERWISE SKIP\n",
    "\n",
    "using CSV;\n",
    "\n",
    "period_length = 10;\n",
    "agent_type_n = 10;\n",
    "time_horizon = 500;\n",
    "\n",
    "type_labels = Vector{String}(undef, agent_type_n + 2)\n",
    "for i in eachindex(type_labels)\n",
    "    type_labels[i] = string(\"f\", i-2)\n",
    "end\n",
    "type_labels[1] = \"beta\"\n",
    "type_labels[2] = \"eta\"\n",
    "\n",
    "model_instance = \"Y$(time_horizon)T$(period_length)N$(agent_type_n)Benchmark\";\n",
    "\n",
    "dogmatism_runs = 1\n",
    "\n",
    "modelrun = NamedTuple[(\n",
    "                degree = 1000.0,\n",
    "                βs = [ (1+.015)^(-10); (1+.001)^(-10)],\n",
    "                ηs = [1.45; 1.0],\n",
    "                weights = [1 0; 0 1])]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Initialise the DICE-like model with cumulative emissions climate module and calibrate to DICE-2016R\n",
    "\n",
    "using Plots\n",
    "using JuMP, Ipopt, BasisMatrices\n",
    "using Parameters\n",
    "\n",
    "T = div(time_horizon, period_length)\n",
    "\n",
    "runnumber = 2\n",
    "solution_file = string(\"run$(runnumber)degree$(Int(10000-10000*modelrun[runnumber].degree))\")\n",
    "\n",
    "print(\"Setting up model run parameters...\")\n",
    "\n",
    "# Heterogeneous types\n",
    "βs = modelrun[runnumber].βs                 # Discount factor\n",
    "ηs = modelrun[runnumber].ηs                 # Elasticity of MU\n",
    "agents = Array{NamedTuple,1}(undef, length(βs))\n",
    "for (i, (β, η)) in enumerate(zip(βs, ηs))\n",
    "   agents[i] = NamedTuple{(:β, :η),Tuple{Float64, Float64}}((β, η))\n",
    "end\n",
    "\n",
    "weights = modelrun[runnumber].weights\n",
    "\n",
    "println(\" done.\")\n",
    "\n",
    "\n",
    "A_seq = zeros(T)\n",
    "\n",
    "print(\"Calculating population path...\")\n",
    "\n",
    "# Population\n",
    "# Taken directly from DICE-2016 GAMS version, September 2016\n",
    "L_seq = ones(T)\n",
    "L_ultimate = 11.500\n",
    "L_seq[1] = 7.403\n",
    "L_adjustment_annual = .02836              # In DICE, L^+/Lbar = (L/Lbar)^(1 - .134), so annual adjustment is 1 - (1-.134)^.2\n",
    "L_adjustment = 1 - (1 - L_adjustment_annual)^period_length\n",
    "for t in 2:T\n",
    "    L_seq[t] = L_seq[t-1] * (L_ultimate / L_seq[t-1])^L_adjustment\n",
    "end\n",
    "\n",
    "println(\" done.\")\n",
    "\n",
    "\n",
    "print(\"Calculating TFP path...\")\n",
    "\n",
    "# TFP sequence\n",
    "# Adapted to match  DICE-2016 GAMS version, September 2016\n",
    "A_init = 5.115\n",
    "A_growth = ones(T)\n",
    "A_growth[1] =  (1 + .0148)^period_length - 1                # Annual growth rate for one year (from DICE's .076 for five years: equal to (1.076)^.2 )\n",
    "A_growth_decay = .005\n",
    "A_seq[1] = A_init\n",
    "\n",
    "for t in 2:T\n",
    "    A_growth[t] = A_growth[1] * exp(-A_growth_decay * period_length * (t - 1))\n",
    "    A_seq[t] = A_seq[t - 1] / (1 - A_growth[t])\n",
    "end\n",
    "\n",
    "println(\" done.\")\n",
    "\n",
    "\n",
    "print(\"Setting up production and utility functions...\")\n",
    "# Gross output\n",
    "α = .3    # Capital share\n",
    "F(A, K, L) = A * K^α * L^(1 - α)\n",
    "dF(A, K, L) = α * F(A, K, L)./K\n",
    "\n",
    "# Labour-augmenting A\n",
    "Ae_seq = A_seq.^(1/(1 - α))\n",
    "\n",
    "# Consumption utility\n",
    "u(c, η) = (η == 1.0) ? log.(c/10) : (c/10).^(1-η)/(1 - η)\n",
    "du(c, η) = (η == 1.0) ? 1.0 ./ ((c/10) * 10) : ((c/10).^(-η))/10\n",
    "\n",
    "println(\" done.\")\n",
    "\n",
    "\n",
    "print(\"Setting up damage functions...\")\n",
    "\n",
    "# Climate damage\n",
    "# Note that the calibration of climate damage is implicitly used to determine the upper bound on the S-grid\n",
    "γ1 = 0.0     # Linear term\n",
    "γ2 = .00236   # Convex term level from DICE, adjusted to GtC\n",
    "γ3 = 2.0     # Convexity\n",
    "\n",
    "# We work with cumulative carbon emissions, so that the preindustrial stock is just zero\n",
    "#    It will then be reasonable to suppose 500 GtC extracted by about 2010, so S_0 = 0.5 TtC\n",
    "S_preind = 0*2.13 # Preindustrial stock of 280 ppm in GtC\n",
    "temperature_sensitivity = 2.0\n",
    "\n",
    "# Damage(S) = 1 - γ1 * S - γ2 * S^γ3\n",
    "Damage(S) = γ1 .* (temperature_sensitivity * S) + γ2 .* (temperature_sensitivity * S).^γ3\n",
    "dDamage(S) = γ1 + γ2 * γ3 .* (temperature_sensitivity * S).^(γ3-1)\n",
    "\n",
    "println(\" done.\")\n",
    "\n",
    "\n",
    "print(\"Calculating carbon intensity path...\")\n",
    "\n",
    "# Carbon intensity\n",
    "E_init = 35.85                           # Note: these emissions at in GtCO2 per year, not GtC!\n",
    "Q_init = 105.5\n",
    "μ_init = .03\n",
    "\n",
    "σ = ones(T)\n",
    "g_σ = ones(T)\n",
    "\n",
    "\n",
    "# Carbon intensity\n",
    "σ[1] = E_init/(Q_init * (1 - μ_init)) * 12 / 44         # The last adjustment makes sure we measure carbon emission intensity in GtC\n",
    "g_σ[1] = -.0152\n",
    "g_σ_decay = -.001\n",
    "    \n",
    "for t = 2:T\n",
    "    g_σ[t] = g_σ[t-1]*(1 + g_σ_decay)^period_length\n",
    "    σ[t] = σ[t-1] * exp(g_σ[t-1]*period_length)\n",
    "end\n",
    "\n",
    "println(\" done.\")\n",
    "\n",
    "\n",
    "print(\"Calculating backstop price path...\")\n",
    "\n",
    "# Backstop price\n",
    "# Taken directly from DICE-2016 GAMS version, September 2016\n",
    "p_backstop_init = 550*44/12                # pB reported in $/tCO2 in Nordhaus. We use $/tC.\n",
    "g_backstop = .005051                       # Backstop price decrease rate. In DICE, .025 per five years. The annual rate equals (1 - .025)^.2\n",
    "p_backstop = ones(T)\n",
    "for t in 1:T\n",
    "    p_backstop[t] = p_backstop_init*(1 - g_backstop)^(period_length * (t - 1))\n",
    "end\n",
    "\n",
    "println(\" done.\")\n",
    "\n",
    "\n",
    "print(\"Setting up abatement cost function...\")\n",
    "        \n",
    "# Abatement cost\n",
    "θ2 = 2.6    # Cost convexity\n",
    "θ1 = ones(T)\n",
    "θ1 = (p_backstop / 1000) .* σ ./ θ2                   # Calibrated so that full decarbonisation implies marginal cost equal to backstop use\n",
    "Abatecost(μ,theta) = theta * μ .^θ2\n",
    "dAbatecost(μ, theta) = theta * θ2 * μ.^(θ2 - 1)\n",
    "\n",
    "println(\" done.\")\n",
    "\n",
    "\n",
    "print(\"Setting up land use emissions...\")\n",
    "E_land_init = 2.6 * 12/44;                 # Land-use change emissions in 2015 in GtC\n",
    "E_land_decay = .115;                       # Rate of decrease of land-use emissions\n",
    "E_land_seq = ones(T);\n",
    "for t in 1:T\n",
    "    E_land_seq[t] = E_land_init*(1 - E_land_decay)^(period_length/5*(t-1));\n",
    "end\n",
    "\n",
    "println(\" done.\")\n",
    "\n",
    "\n",
    "print(\"Setting up net output and state evolution functions...\")\n",
    "\n",
    "# Net output\n",
    "# Taken directly from DICE-2016 GAMS version, September 2016\n",
    "Y(A, K, L, μ, S, θ) = (1 - Damage(S) - Abatecost(μ,θ)) .* F(A, K, L)\n",
    "dY(A, K, L, μ, S, θ) = (1 - Damage(S) - Abatecost(μ,θ)) .* dF(A, K, L)\n",
    "\n",
    "# Capital and carbon dynamics\n",
    "δ = 1 - (1 - .10)^period_length    # Capital depreciation using DICE framework, converted into a depreciation rate for our time period\n",
    "K_next(K, Y, C) = (1 - δ) * K + period_length*(Y - C)\n",
    "\n",
    "# These function convert capital stocks into effective capital stocks (normalised by TFP, population) and back\n",
    "eff(X, t) = (X / (L_seq[t] * Ae_seq[t]))\n",
    "uneff(X, t) = X * (L_seq[t] * Ae_seq[t])\n",
    "\n",
    "S_next(S, F, μ, σ) = S + period_length * ((σ * (1 - μ) * F) / 1000 + E_land_seq[t] / 1000)\n",
    "\n",
    "println(\" done.\")\n",
    "\n",
    "\n",
    "print(\"Setting up the approximation grid...\")\n",
    "\n",
    "# Set the endpoints of approximation interval\n",
    "# Note that the capital stock here is effective capital stock\n",
    "Ke_max = ceil(1.1*(α * period_length / δ)^(1/(1 - α)))           # 10 percent above maximal Golden Rule capital\n",
    "S_max = ceil(((.2/γ2)^.5)/temperature_sensitivity)              # .2 refers to a level of climate damages we expect to never reach\n",
    "\n",
    "# In practice the simulations require the capital stock upper cap to be high enough to work.\n",
    "# Below are some values which have worked.  In case of numerical instability, first thing to\n",
    "# try is to increase the capital stock upper bound.\n",
    "\n",
    "# These values work for run 1\n",
    "#a̲ = [Ke_max/20 .450]                            # left endpoint\n",
    "#a̅= [1.8*Ke_max S_max]                            # right endpoint\n",
    "\n",
    "# These values work for run 2\n",
    "a̲ = [Ke_max/20 .450]                            # left endpoint\n",
    "a̅= [5*Ke_max S_max]                            # right endpoint\n",
    "\n",
    "# This has worked well.  Lower orders have sometimes led to numerical instability.\n",
    "n = [24 10]                             # order of approximation\n",
    "\n",
    "# Set up basis matrices for the value function approximations\n",
    "basis = Basis(ChebParams(n[1], a̲[1], a̅[1]), ChebParams(n[2], a̲[2], a̅[2]))    # define basis\n",
    "stategrid, (Kegrid, Sgrid) = nodes(basis)\n",
    "\n",
    "Ψ = BasisMatrix(basis, Expanded(), stategrid).vals[1]\n",
    "Ψ_S = BasisMatrix(basis, Expanded(), stategrid, [0 1]).vals[1]\n",
    "Ψ_K = BasisMatrix(basis, Expanded(), stategrid, [1 0]).vals[1]\n",
    "\n",
    "println(\" done.\")\n",
    "\n",
    "\n",
    "\n",
    "print(\"Setting up approximation tables...\")\n",
    "\n",
    "C = zeros(T, length(βs))\n",
    "V = Array{Function, 2}(undef, T, length(βs))\n",
    "Φ = Array{Array{Float64}}(undef, T, length(βs))\n",
    "optimal_mu = Array{Array{Float64}}(undef, T, length(βs))\n",
    "optimal_cons = similar(optimal_mu)\n",
    "Values = Array{Array{Float64}}(undef, T, length(βs))\n",
    "\n",
    "println(\" done.\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up the parameter vectors and functions used to solve the model\n",
    "\n",
    "DICE_parameters = (A = A_seq, Ae = Ae_seq, E_land = E_land_seq, L_seq = L_seq, ηs = ηs, βs = βs, weights = weights, θ1 = θ1, δ = δ, σ = σ, T = T)\n",
    "approx_parameters = (a̲ = a̲, a̅ = a̅, n = n, period_length = period_length)\n",
    "algo_parameters = (tries = 30, initial_control = [100.0 0.5], verbosity = 0)\n",
    "debug_parameters = (tries = 1, initial_control = [100.0 0.5], verbosity = 5)\n",
    "functions = ( u = u, du = du, F = F, dF = dF, Damage = Damage, dDamage = dDamage, Abatecost = Abatecost, dAbatecost = dAbatecost, K_next = K_next, S_next = S_next, eff = eff, uneff = uneff )\n",
    "\n",
    "function solve_coefficients(agent, time, Φ, stategrid, V, opt_cons, opt_mu, objective_weights, DICE_parameters, approx_parameters, functions, verbosity)\n",
    "\n",
    "    # Given a set of values computed at the approximation nodes, this function computes the coefficients for a Chebyshev approximation.\n",
    "    # NB. This is not used in the current implementation, rather the coefficients are obtained using a matrix inversion later.\n",
    "    \n",
    "    time_1 = time_ns();\n",
    "    print_debug = false\n",
    "    @unpack A, L_seq, E_land, ηs, βs, weights, θ1, δ, σ, T, Ae = DICE_parameters\n",
    "    @unpack a̲, a̅, n = approx_parameters\n",
    "    @unpack u, du, F, dF, Damage, dDamage, Abatecost, dAbatecost, K_next, S_next, eff, uneff = functions\n",
    "\n",
    "    chebpols = zeros(size(stategrid, 1), n[1]*n[2])\n",
    "    chebpols_K = similar(chebpols)\n",
    "    chebpols_S = similar(chebpols)\n",
    "    chebpols_KK = similar(chebpols)\n",
    "    chebpols_next = zeros(length(βs), size(stategrid, 1), n[1]*n[2])\n",
    "    chebpols_next_K = similar(chebpols_next)\n",
    "    chebpols_next_S = similar(chebpols_next)\n",
    "    \n",
    "    for j = 1:size(stategrid, 1)\n",
    "        chebpols[j, :] = (BasisMatrix(Basis(ChebParams(n[1], a̲[1], a̅[1]), ChebParams(n[2], a̲[2], a̅[2])), \n",
    "            Expanded(), [stategrid[j,1], stategrid[j,2]]).vals[1])\n",
    "        chebpols_K[j, :] = (BasisMatrix(Basis(ChebParams(n[1], a̲[1], a̅[1]), ChebParams(n[2], a̲[2], a̅[2])), \n",
    "            Expanded(), [stategrid[j,1], stategrid[j,2]], [1 0]).vals[1])\n",
    "        chebpols_KK[j, :] = (BasisMatrix(Basis(ChebParams(n[1], a̲[1], a̅[1]), ChebParams(n[2], a̲[2], a̅[2])), \n",
    "            Expanded(), [stategrid[j,1], stategrid[j,2]], [2 0]).vals[1])\n",
    "        chebpols_S[j, :] = (BasisMatrix(Basis(ChebParams(n[1], a̲[1], a̅[1]), ChebParams(n[2], a̲[2], a̅[2])), \n",
    "            Expanded(), [stategrid[j,1], stategrid[j,2]], [0 1]).vals[1])\n",
    "        if (time < T)\n",
    "            for m = 1:length(βs)\n",
    "                chebpols_next[m, j, :] = (BasisMatrix(Basis(ChebParams(n[1], a̲[1], a̅[1]), ChebParams(n[2], a̲[2], a̅[2])), \n",
    "                    Expanded(), [eff.((1 - δ) * uneff(stategrid[j,1], time) + period_length*(Y(A[time], uneff(stategrid[j,1], time), L_seq[time], opt_mu[j], stategrid[j,2], θ1[time]) - opt_cons[j]), time + 1), stategrid[j,2] + period_length*(σ[time] * (1 - opt_mu[j]) * F(A[time], uneff(stategrid[j,1], time), L_seq[time]) + E_land[time])/1000]).vals[1])\n",
    "                chebpols_next_K[m, j, :] = (BasisMatrix(Basis(ChebParams(n[1], a̲[1], a̅[1]), ChebParams(n[2], a̲[2], a̅[2])), \n",
    "                    Expanded(), [eff.((1 - δ) * uneff(stategrid[j,1], time) + period_length*(Y(A[time], uneff(stategrid[j,1], time), L_seq[time], opt_mu[j], stategrid[j,2], θ1[time]) - opt_cons[j]), time + 1), stategrid[j,2] + period_length*(σ[time] * (1 - opt_mu[j]) * F(A[time], uneff(stategrid[j,1], time), L_seq[time]) + E_land[time])/1000], [1 0]).vals[1])\n",
    "                chebpols_next_S[m, j, :] = (BasisMatrix(Basis(ChebParams(n[1], a̲[1], a̅[1]), ChebParams(n[2], a̲[2], a̅[2])), \n",
    "                    Expanded(), [eff.((1 - δ) * uneff(stategrid[j,1], time) + period_length*(Y(A[time], uneff(stategrid[j,1], time), L_seq[time], opt_mu[j], stategrid[j,2], θ1[time]) - opt_cons[j]), time + 1), stategrid[j,2] + period_length*(σ[time] * (1 - opt_mu[j]) * F(A[time], uneff(stategrid[j,1], time), L_seq[time]) + E_land[time])/1000], [0 1]).vals[1])\n",
    "            end\n",
    "        end\n",
    "    end\n",
    "\n",
    "    model = Model(with_optimizer(Ipopt.Optimizer, max_iter=20000, print_level=5, bound_relax_factor=0.0, tol=1e-6, ))\n",
    "\n",
    "    @NLparameter(model, A == A[time])\n",
    "    @NLparameter(model, E_land == E_land[time])\n",
    "    @NLparameter(model, L == L_seq[time])\n",
    "    @NLparameter(model, σ == σ[time])\n",
    "    @NLparameter(model, η == ηs[agent])\n",
    "    @NLparameter(model, β == βs[agent])\n",
    "    @NLparameter(model, Ke[h = 1:size(stategrid,1)] == stategrid[h,1])\n",
    "    @NLparameter(model, S[h = 1:size(stategrid,1)] == stategrid[h,2])\n",
    "    @NLparameter(model, w[h = 1:length(βs)] == weights[agent, h])\n",
    "    @NLparameter(model, eff == Ae[time] * L_seq[time])\n",
    "    if (time < T)\n",
    "        @NLparameter(model, phi[h = 1:length(βs), k = 1:n[1]*n[2]] == Φ[time+1, h][k])\n",
    "        @NLparameter(model, eff_next == Ae[time + 1] * L_seq[time + 1])\n",
    "    end\n",
    "    @NLparameter(model, θ1 == θ1[time])\n",
    "    @NLparameter(model, opt_cons[h = 1:length(opt_cons)] == opt_cons[h])\n",
    "    @NLparameter(model, opt_mu[h = 1:length(opt_mu)] == opt_mu[h])\n",
    "\n",
    "    Ψ = BasisMatrix(basis, Expanded(), stategrid).vals[1]\n",
    "    temp_coefficients = Ψ \\ V\n",
    "    \n",
    "    @variable(model, phi_new[k = 1:n[1]*n[2]], start = temp_coefficients[k])\n",
    "\n",
    "    JuMP.register(model, :F, 3, F, autodiff=true)\n",
    "    JuMP.register(model, :Damage, 1, Damage, autodiff=true)\n",
    "    JuMP.register(model, :Abatecost, 2, Abatecost, autodiff=true)\n",
    "    JuMP.register(model, :Y, 6, Y, autodiff=true)\n",
    "    JuMP.register(model, :du, 2, du, autodiff=true)\n",
    "    JuMP.register(model, :dF, 3, dF, autodiff=true)\n",
    "    JuMP.register(model, :dDamage, 1, dDamage, autodiff=true)\n",
    "    JuMP.register(model, :dAbatecost, 2, dAbatecost, autodiff=true)\n",
    "\n",
    "    @NLparameter(model, chebpols[j in 1:size(stategrid, 1), h in 1:n[1]*n[2]] == chebpols[j, h])\n",
    "    @NLparameter(model, chebpols_K[j in 1:size(stategrid, 1), h in 1:n[1]*n[2]] == chebpols_K[j, h])\n",
    "    @NLparameter(model, chebpols_S[j in 1:size(stategrid, 1), h in 1:n[1]*n[2]] == chebpols_S[j, h])\n",
    "    @NLparameter(model, chebpols_KK[j in 1:size(stategrid, 1), h in 1:n[1]*n[2]] == chebpols_KK[j, h])\n",
    "    #@NLparameter(model, chebpols_next[i in 1:length(βs), j in 1:size(stategrid, 1), h in 1:n[1]*n[2]] == chebpols_next[i, j, h])\n",
    "    @NLparameter(model, chebpols_next_K[i in 1:length(βs), j in 1:size(stategrid, 1), h in 1:n[1]*n[2]] == chebpols_next_K[i, j, h])\n",
    "    @NLparameter(model, chebpols_next_S[i in 1:length(βs), j in 1:size(stategrid, 1), h in 1:n[1]*n[2]] == chebpols_next_S[i, j, h])\n",
    "    \n",
    "    @variable(model, V_app[j in 1:size(stategrid, 1)], start = V[j])\n",
    "    @variable(model, 0.0 <= V_K_app[j in 1:size(stategrid, 1)])\n",
    "    @variable(model, V_S_app[j in 1:size(stategrid, 1)] <= 0.0)\n",
    "    @variable(model, V_KK_app[j in 1:size(stategrid, 1)] <= 0.0)\n",
    "    @NLconstraint(model, [j = 1:size(stategrid, 1)], V_app[j] == sum(chebpols[j, h] * phi_new[h] for h in 1:n[1]*n[2]))\n",
    "    @NLconstraint(model, [j = 1:size(stategrid, 1)], V_K_app[j] == sum(chebpols_K[j, h] * phi_new[h] for h in 1:n[1]*n[2]))\n",
    "    @NLconstraint(model, [j = 1:size(stategrid, 1)], V_S_app[j] == sum(chebpols_S[j, h] * phi_new[h] for h in 1:n[1]*n[2]))\n",
    "    @NLconstraint(model, [j = 1:size(stategrid, 1)], V_KK_app[j] == sum(chebpols_KK[j, h] * phi_new[h] for h in 1:n[1]*n[2]))\n",
    "    \n",
    "    if time < T\n",
    "        #@variable(model, V_next_app[m = 1:length(βs), j in 1:size(stategrid, 1)])\n",
    "        @variable(model, V_K_next_app[m = 1:length(βs), j in 1:size(stategrid, 1)])\n",
    "        @variable(model, V_S_next_app[m = 1:length(βs), j in 1:size(stategrid, 1)])\n",
    "        #@NLconstraint(model, [m = 1:length(βs), j = 1:size(stategrid, 1)], V_next_app[m,j] == sum(chebpols_next[m, j, h] * phi[m, h] for h in 1:n[1]*n[2]))\n",
    "        @NLconstraint(model, [m = 1:length(βs), j = 1:size(stategrid, 1)], V_K_next_app[m,j] == sum(chebpols_next_K[m, j, h] * phi[m, h] for h in 1:n[1]*n[2]))\n",
    "        @NLconstraint(model, [m = 1:length(βs), j = 1:size(stategrid, 1)], V_S_next_app[m,j] == sum(chebpols_next_S[m, j, h] * phi[m, h] for h in 1:n[1]*n[2]))\n",
    "    end\n",
    "    \n",
    "    if print_debug\n",
    "        println(\"JuMP model set up, except for objective.\")\n",
    "    end\n",
    "    #time_2 = time_ns()\n",
    "    #println(\"$(time_2 - time_1)\")\n",
    "    #time_1 = time_ns()\n",
    "    \n",
    "    # Minimise by choosing the period t coefficients:\n",
    "    #   summing over all the grid points\n",
    "    #       Sq deviation of Value less Approximate value \n",
    "    #           Here, calculation Approximation by V_current using the coefficients\n",
    "    #       Sq deviation of V_K today, by env theorem only dep. on less expected next period V_K and V_S\n",
    "    #           Summing over future types\n",
    "    #           Calculating V_K using future coefficients\n",
    "    #       Sq deviation of V_S today, similarly using env theorem\n",
    "    \n",
    "    if (time < T)\n",
    "        if print_debug\n",
    "            println(\"Optimising for t < T.\")\n",
    "        end\n",
    "        @NLobjective(model, Min, sum(\n",
    "                objective_weights[1] * (V[j] - V_app[j])^2               # residual between optimised value and the new value fn\n",
    "                + objective_weights[2] * ( V_K_app[j] # V_K_now_aux\n",
    "                    - β * eff * sum(w[m] * ((1 - δ + period_length * (1 - Damage(S[j]) - Abatecost(opt_mu[j], θ1)) * dF(A, eff * Ke[j], L) ) * V_K_next_app[m, j] / eff_next # V_K_next\n",
    "                        + period_length * σ * (1 - opt_mu[j]) * dF(A, eff * Ke[j], L) / 1000 * V_S_next_app[m, j] # V_S_next\n",
    "                        ) for m = 1:length(βs)) # end of sum\n",
    "                    )^2\n",
    "                + objective_weights[3] * ( V_S_app[j] # V_S_now\n",
    "                    - β * sum(w[m] * ((-period_length * dDamage(S[j]) * F(A, eff * Ke[j], L)) * V_K_next_app[m,j] / eff_next + # V_K_next\n",
    "                          V_S_next_app[m,j] # V_S_next\n",
    "                        ) for m = 1:length(βs))\n",
    "                    )^2\n",
    "            for j = 1:length(Ke)) # closes the sum over j\n",
    "         ) # closes the objective call\n",
    "        else # In the last period there is no continuation, just the current capital and output will be consumed\n",
    "        @NLobjective(model, Min, sum( # THESE NEED CHANGING TOO FOR PROD FN SPEC CHANGE\n",
    "                objective_weights[1]*(V[j] - V_app[j])^2            # residual between optimised value and the new value fn\n",
    "                + objective_weights[2] * ( V_K_app[j] # Residual between current V_K...\n",
    "                    - eff * du(opt_cons[j]/L, η) * ((1 - δ + (1 - Damage(S[j]) - Abatecost(opt_mu[j], θ1)) * dF(A, eff * Ke[j], L) )))^2\n",
    "                + objective_weights[3] * ( V_S_app[j] # V_S_now_aux\n",
    "                    - du(opt_cons[j]/L, η) * (-dDamage(S[j]) * F(A, eff * Ke[j], L)))^2\n",
    "            for j = 1:length(Ke)) # closes the sum over j\n",
    "         ) # closes the objective call        \n",
    "    end\n",
    "    if print_debug\n",
    "        println(\"JuMP model objective set up.\")\n",
    "    end\n",
    "    #time_2 = time_ns()\n",
    "    #println(\"$(time_2 - time_1)\")\n",
    "    #time_1 = time_ns()\n",
    "    \n",
    "    flush(stdout)\n",
    "    optimize!(model);\n",
    "    if print_debug\n",
    "        println(\"Model solved.\")\n",
    "    end\n",
    "    #time_2 = time_ns()\n",
    "    #println(\"$(time_2 - time_1)\")\n",
    "\n",
    "    return objective_value(model), value.(phi_new), termination_status(model)\n",
    "    \n",
    "end\n",
    "\n",
    "\n",
    "function solve_decision(state, agent, time, Φ, DICE_parameters, approx_parameters, initial_control, functions, tolerance, verbosity)\n",
    "    \n",
    "    # Given a state and a next period value function, this function computes the optimal decision\n",
    "\n",
    "    print_debug = false\n",
    "    Ke, S = state\n",
    "    @unpack A, L_seq, ηs, βs, weights, θ1, δ, σ, T, Ae, E_land = DICE_parameters\n",
    "    @unpack a̲, a̅, n = approx_parameters\n",
    "    @unpack u, du, F, dF, Damage, dDamage, Abatecost, dAbatecost, K_next, S_next = functions\n",
    "    \n",
    "    if print_debug\n",
    "        println(\"Parameters unpacked.\")\n",
    "    end\n",
    "        \n",
    "    #model = Model(solver = IpoptSolver(max_iter=10000, print_level=0, bound_relax_factor=0.0, tol=1e-6))\n",
    "    model = Model(with_optimizer(Ipopt.Optimizer, max_iter=10000, print_level=verbosity, bound_relax_factor=0.0, tol=tolerance, acceptable_constr_viol_tol=1e-5, acceptable_tol=1e-8, acceptable_iter=30))\n",
    "    if print_debug\n",
    "        println(\"Model set up.\")\n",
    "    end\n",
    "\n",
    "    @NLparameter(model, A == A[time])\n",
    "    @NLparameter(model, E_land == E_land_seq[time])\n",
    "    @NLparameter(model, L == L_seq[time])\n",
    "    @NLparameter(model, σ == σ[time])\n",
    "    @NLparameter(model, η == ηs[agent])\n",
    "    @NLparameter(model, β == βs[agent])\n",
    "    @NLparameter(model, Ke == Ke)\n",
    "    @NLparameter(model, S == S)\n",
    "    @NLparameter(model, w[h = 1:length(βs)] == weights[agent, h])\n",
    "    @NLparameter(model, phi[h = 1:length(βs), k = 1:n[1]*n[2]] == Φ[time+1, h][k])\n",
    "    @NLparameter(model, θ1 == θ1[time])\n",
    "    @NLparameter(model, eff == Ae[time] * L_seq[time])\n",
    "    @NLparameter(model, eff_next == Ae[time + 1] * L_seq[time + 1])\n",
    "    if print_debug\n",
    "        println(\"Parameters loaded into model.\")\n",
    "    end\n",
    "\n",
    "    @variable(model, cons >= 0.00001, start=initial_control[1])\n",
    "    #@variable(model, aux_K_next => 0, start = K)\n",
    "    #@variable(model, 0.0 <= cons <= , start=10)\n",
    "    \n",
    "    \n",
    "    @variable(model, 1e-4 <= mu <= 1, start=initial_control[2])\n",
    "    if print_debug\n",
    "        println(\"Variables loaded into model.\")\n",
    "    end\n",
    "\n",
    "    JuMP.register(model, :u, 2, u, autodiff=true)\n",
    "    JuMP.register(model, :F, 3, F, autodiff=true)\n",
    "    JuMP.register(model, :Damage, 1, Damage, autodiff=true)\n",
    "    JuMP.register(model, :Abatecost, 2, Abatecost, autodiff=true)\n",
    "    JuMP.register(model, :Y, 6, Y, autodiff=true)\n",
    "    if print_debug\n",
    "        println(\"Functions defined in model.\")\n",
    "    end\n",
    "    \n",
    "    @NLexpression(model, Ke_next, ((1 - δ) * eff * Ke + period_length * (Y(A, eff * Ke, L, mu, S, θ1) - cons)) / eff_next);\n",
    "    @NLexpression(model, S_next, S + period_length * (σ * (1 - mu) * F(A, eff * Ke, L) + E_land) / 1000);\n",
    "    @NLconstraint(model, Ke_positive, a̲[1] + 1e-6 <= Ke_next <= a̅[1])\n",
    "    if print_debug\n",
    "        println(\"Constraints loaded into model.\")\n",
    "    end\n",
    "\n",
    "    \n",
    "    function V_next(Ke_next, S_next, phi...)\n",
    "        (BasisMatrix(Basis(ChebParams(n[1], a̲[1], a̅[1]), ChebParams(n[2], a̲[2], a̅[2])), \n",
    "            Expanded(), [Ke_next, S_next]).vals[1] * collect(phi))[1]\n",
    "        #(BasisMatrix(Basis(SplineParams(Kgrid0, 0, 3), SplineParams(Sgrid0, 0, 3)), \n",
    "        #Expanded(), [K_next, S_next]).vals[1] * collect(phi))[1]\n",
    "    end\n",
    "    \n",
    "    JuMP.register(model, :V_next, n[1]*n[2]+2, V_next, autodiff=true)\n",
    "    if print_debug\n",
    "        println(\"V_next defined in model.\")\n",
    "    end\n",
    "    \n",
    "    # Greetings if you have made it this far! Below you see a downside of using JuMP.\n",
    "    @NLobjective(model, Max, u(cons/L, η)*L + β * sum(w[m] * V_next(Ke_next, S_next,\n",
    "                    phi[m, 1], phi[m, 2], phi[m, 3], phi[m, 4], phi[m, 5], phi[m, 6], phi[m, 7], phi[m, 8], phi[m, 9], phi[m, 10],\n",
    "                    phi[m, 11], phi[m, 12], phi[m, 13], phi[m, 14], phi[m, 15], phi[m, 16], phi[m, 17], phi[m, 18], phi[m, 19], phi[m, 20],\n",
    "                    phi[m, 21], phi[m, 22], phi[m, 23], phi[m, 24], phi[m, 25], phi[m, 26], phi[m, 27], phi[m, 28], phi[m, 29], phi[m, 30],\n",
    "                    phi[m, 31], phi[m, 32], phi[m, 33], phi[m, 34], phi[m, 35], phi[m, 36], phi[m, 37], phi[m, 38], phi[m, 39], phi[m, 40],\n",
    "                    phi[m, 41], phi[m, 42], phi[m, 43], phi[m, 44], phi[m, 45], phi[m, 46], phi[m, 47], phi[m, 48], phi[m, 49], phi[m, 50],\n",
    "                    phi[m, 51], phi[m, 52], phi[m, 53], phi[m, 54], phi[m, 55], phi[m, 56], phi[m, 57], phi[m, 58], phi[m, 59], phi[m, 60],\n",
    "                    phi[m, 61], phi[m, 62], phi[m, 63], phi[m, 64], phi[m, 65], phi[m, 66], phi[m, 67], phi[m, 68], phi[m, 69], phi[m, 70],\n",
    "                    phi[m, 71], phi[m, 72], phi[m, 73], phi[m, 74], phi[m, 75], phi[m, 76], phi[m, 77], phi[m, 78], phi[m, 79], phi[m, 80],\n",
    "                    phi[m, 81], phi[m, 82], phi[m, 83], phi[m, 84], phi[m, 85], phi[m, 86], phi[m, 87], phi[m, 88], phi[m, 89], phi[m, 90],\n",
    "                    phi[m, 91], phi[m, 92], phi[m, 93], phi[m, 94], phi[m, 95], phi[m, 96], phi[m, 97], phi[m, 98], phi[m, 99], phi[m, 100],\n",
    "                    phi[m, 101], phi[m, 102], phi[m, 103], phi[m, 104], phi[m, 105], phi[m, 106], phi[m, 107], phi[m, 108], phi[m, 109], phi[m, 110],\n",
    "                    phi[m, 111], phi[m, 112], phi[m, 113], phi[m, 114], phi[m, 115], phi[m, 116], phi[m, 117], phi[m, 118], phi[m, 119], phi[m, 120],\n",
    "                    phi[m, 121], phi[m, 122], phi[m, 123], phi[m, 124], phi[m, 125], phi[m, 126], phi[m, 127], phi[m, 128], phi[m, 129], phi[m, 130],\n",
    "                    phi[m, 131], phi[m, 132], phi[m, 133], phi[m, 134], phi[m, 135], phi[m, 136], phi[m, 137], phi[m, 138], phi[m, 139], phi[m, 140],\n",
    "                    phi[m, 141], phi[m, 142], phi[m, 143], phi[m, 144], phi[m, 145], phi[m, 146], phi[m, 147], phi[m, 148], phi[m, 149], phi[m, 150],\n",
    "                    phi[m, 151], phi[m, 152], phi[m, 153], phi[m, 154], phi[m, 155], phi[m, 156], phi[m, 157], phi[m, 158], phi[m, 159], phi[m, 160],\n",
    "                    phi[m, 161], phi[m, 162], phi[m, 163], phi[m, 164], phi[m, 165], phi[m, 166], phi[m, 167], phi[m, 168], phi[m, 169], phi[m, 170],\n",
    "                    phi[m, 171], phi[m, 172], phi[m, 173], phi[m, 174], phi[m, 175], phi[m, 176], phi[m, 177], phi[m, 178], phi[m, 179], phi[m, 180],\n",
    "                    phi[m, 181], phi[m, 182], phi[m, 183], phi[m, 184], phi[m, 185], phi[m, 186], phi[m, 187], phi[m, 188], phi[m, 189], phi[m, 190],\n",
    "                    phi[m, 191], phi[m, 192], phi[m, 193], phi[m, 194], phi[m, 195], phi[m, 196], phi[m, 197], phi[m, 198], phi[m, 199], phi[m, 200],\n",
    "                    phi[m, 201], phi[m, 202], phi[m, 203], phi[m, 204], phi[m, 205], phi[m, 206], phi[m, 207], phi[m, 208], phi[m, 209], phi[m, 210],\n",
    "                    phi[m, 211], phi[m, 212], phi[m, 213], phi[m, 214], phi[m, 215], phi[m, 216], phi[m, 217], phi[m, 218], phi[m, 219], phi[m, 220],\n",
    "                    phi[m, 221], phi[m, 222], phi[m, 223], phi[m, 224], phi[m, 225], phi[m, 226], phi[m, 227], phi[m, 228], phi[m, 229], phi[m, 230],\n",
    "                    phi[m, 231], phi[m, 232], phi[m, 233], phi[m, 234], phi[m, 235], phi[m, 236], phi[m, 237], phi[m, 238], phi[m, 239], phi[m, 240]#=,\n",
    "                    phi[m, 241], phi[m, 242], phi[m, 243], phi[m, 244], phi[m, 245], phi[m, 246], phi[m, 247], phi[m, 248], phi[m, 249], phi[m, 250],\n",
    "                    phi[m, 251], phi[m, 252], phi[m, 253], phi[m, 254], phi[m, 255], phi[m, 256], phi[m, 257], phi[m, 258], phi[m, 259], phi[m, 260],\n",
    "                    phi[m, 261], phi[m, 262], phi[m, 263], phi[m, 264], phi[m, 265], phi[m, 266], phi[m, 267], phi[m, 268], phi[m, 269], phi[m, 270],\n",
    "                    phi[m, 271], phi[m, 272], phi[m, 273], phi[m, 274], phi[m, 275], phi[m, 276], phi[m, 277], phi[m, 278], phi[m, 279], phi[m, 280],\n",
    "                    phi[m, 281], phi[m, 282], phi[m, 283], phi[m, 284], phi[m, 285], phi[m, 286], phi[m, 287], phi[m, 288], phi[m, 289], phi[m, 290],\n",
    "                    phi[m, 291], phi[m, 292], phi[m, 293], phi[m, 294], phi[m, 295], phi[m, 296], phi[m, 297], phi[m, 298], phi[m, 299], phi[m, 300],\n",
    "                    phi[m, 301], phi[m, 302], phi[m, 303], phi[m, 304], phi[m, 305], phi[m, 306], phi[m, 307], phi[m, 308], phi[m, 309], phi[m, 310],\n",
    "                    phi[m, 311], phi[m, 312], phi[m, 313], phi[m, 314], phi[m, 315], phi[m, 316], phi[m, 317], phi[m, 318], phi[m, 319], phi[m, 320],\n",
    "                    phi[m, 321], phi[m, 322], phi[m, 323], phi[m, 324], phi[m, 325], phi[m, 326], phi[m, 327], phi[m, 328], phi[m, 329], phi[m, 330],\n",
    "                    phi[m, 331], phi[m, 332], phi[m, 333], phi[m, 334], phi[m, 335], phi[m, 336], phi[m, 337], phi[m, 338], phi[m, 339], phi[m, 340],\n",
    "                    phi[m, 341], phi[m, 342], phi[m, 343], phi[m, 344], phi[m, 345], phi[m, 346], phi[m, 347], phi[m, 348], phi[m, 349], phi[m, 350],\n",
    "                    phi[m, 351], phi[m, 352], phi[m, 353], phi[m, 354], phi[m, 355], phi[m, 356], phi[m, 357], phi[m, 358], phi[m, 359], phi[m, 360],\n",
    "                    phi[m, 361], phi[m, 362], phi[m, 363], phi[m, 364], phi[m, 365], phi[m, 366], phi[m, 367], phi[m, 368], phi[m, 369], phi[m, 370],\n",
    "                    phi[m, 371], phi[m, 372], phi[m, 373], phi[m, 374], phi[m, 375], phi[m, 376], phi[m, 377], phi[m, 378], phi[m, 379], phi[m, 380],\n",
    "                    phi[m, 381], phi[m, 382], phi[m, 383], phi[m, 384], phi[m, 385], phi[m, 386], phi[m, 387], phi[m, 388], phi[m, 389], phi[m, 390],\n",
    "                    phi[m, 391], phi[m, 392], phi[m, 393], phi[m, 394], phi[m, 395], phi[m, 396], phi[m, 397], phi[m, 398], phi[m, 399], phi[m, 400],\n",
    "                    phi[m, 401], phi[m, 402], phi[m, 403], phi[m, 404], phi[m, 405], phi[m, 406], phi[m, 407], phi[m, 408], phi[m, 409], phi[m, 410],\n",
    "                    phi[m, 411], phi[m, 412], phi[m, 413], phi[m, 414], phi[m, 415], phi[m, 416], phi[m, 417], phi[m, 418], phi[m, 419], phi[m, 420],\n",
    "                    phi[m, 421], phi[m, 422], phi[m, 423], phi[m, 424], phi[m, 425], phi[m, 426], phi[m, 427], phi[m, 428], phi[m, 429], phi[m, 430],\n",
    "                    phi[m, 431], phi[m, 432], phi[m, 433], phi[m, 434], phi[m, 435], phi[m, 436], phi[m, 437], phi[m, 438], phi[m, 439], phi[m, 440],\n",
    "                    phi[m, 441], phi[m, 442], phi[m, 443], phi[m, 444], phi[m, 445], phi[m, 446], phi[m, 447], phi[m, 448], phi[m, 449], phi[m, 450],\n",
    "                    phi[m, 451], phi[m, 452], phi[m, 453], phi[m, 454], phi[m, 455], phi[m, 456], phi[m, 457], phi[m, 458], phi[m, 459], phi[m, 460],\n",
    "                    phi[m, 461], phi[m, 462], phi[m, 463], phi[m, 464], phi[m, 465], phi[m, 466], phi[m, 467], phi[m, 468], phi[m, 469], phi[m, 470],\n",
    "                    phi[m, 471], phi[m, 472], phi[m, 473], phi[m, 474], phi[m, 475], phi[m, 476], phi[m, 477], phi[m, 478], phi[m, 479], phi[m, 480],\n",
    "                    phi[m, 481], phi[m, 482], phi[m, 483], phi[m, 484], phi[m, 485], phi[m, 486], phi[m, 487], phi[m, 488], phi[m, 489], phi[m, 490],\n",
    "                    phi[m, 491], phi[m, 492], phi[m, 493], phi[m, 494], phi[m, 495], phi[m, 496], phi[m, 497], phi[m, 498], phi[m, 499], phi[m, 500],\n",
    "                    phi[m, 501], phi[m, 502], phi[m, 503], phi[m, 504], phi[m, 505], phi[m, 506], phi[m, 507], phi[m, 508], phi[m, 509], phi[m, 510],\n",
    "                    phi[m, 511], phi[m, 512], phi[m, 513], phi[m, 514], phi[m, 515], phi[m, 516], phi[m, 517], phi[m, 518], phi[m, 519], phi[m, 520],\n",
    "                    phi[m, 521], phi[m, 522], phi[m, 523], phi[m, 524], phi[m, 525], phi[m, 526], phi[m, 527], phi[m, 528], phi[m, 529], phi[m, 530],\n",
    "                    phi[m, 531], phi[m, 532], phi[m, 533], phi[m, 534], phi[m, 535], phi[m, 536], phi[m, 537], phi[m, 538], phi[m, 539], phi[m, 540],\n",
    "                    phi[m, 541], phi[m, 542], phi[m, 543], phi[m, 544], phi[m, 545], phi[m, 546], phi[m, 547], phi[m, 548], phi[m, 549], phi[m, 550],\n",
    "                    phi[m, 551], phi[m, 552], phi[m, 553], phi[m, 554], phi[m, 555], phi[m, 556], phi[m, 557], phi[m, 558], phi[m, 559], phi[m, 560],\n",
    "                    phi[m, 561], phi[m, 562], phi[m, 563], phi[m, 564], phi[m, 565], phi[m, 566], phi[m, 567], phi[m, 568], phi[m, 569], phi[m, 570],\n",
    "                    phi[m, 571], phi[m, 572], phi[m, 573], phi[m, 574], phi[m, 575], phi[m, 576], phi[m, 577], phi[m, 578], phi[m, 579], phi[m, 580],\n",
    "                    phi[m, 581], phi[m, 582], phi[m, 583], phi[m, 584], phi[m, 585], phi[m, 586], phi[m, 587], phi[m, 588], phi[m, 589], phi[m, 590],\n",
    "                    phi[m, 591], phi[m, 592], phi[m, 593], phi[m, 594], phi[m, 595], phi[m, 596], phi[m, 597], phi[m, 598], phi[m, 599], phi[m, 600],\n",
    "                    phi[m, 601], phi[m, 602], phi[m, 603], phi[m, 604], phi[m, 605], phi[m, 606], phi[m, 607], phi[m, 608], phi[m, 609], phi[m, 610],\n",
    "                    phi[m, 611], phi[m, 612], phi[m, 613], phi[m, 614], phi[m, 615], phi[m, 616], phi[m, 617], phi[m, 618], phi[m, 619], phi[m, 620],\n",
    "                    phi[m, 621], phi[m, 622], phi[m, 623], phi[m, 624], phi[m, 625], phi[m, 626], phi[m, 627], phi[m, 628], phi[m, 629], phi[m, 630],\n",
    "                    phi[m, 631], phi[m, 632], phi[m, 633], phi[m, 634], phi[m, 635], phi[m, 636], phi[m, 637], phi[m, 638], phi[m, 639], phi[m, 640],\n",
    "                    phi[m, 641], phi[m, 642], phi[m, 643], phi[m, 644], phi[m, 645], phi[m, 646], phi[m, 647], phi[m, 648], phi[m, 649], phi[m, 650],\n",
    "                    phi[m, 651], phi[m, 652], phi[m, 653], phi[m, 654], phi[m, 655], phi[m, 656], phi[m, 657], phi[m, 658], phi[m, 659], phi[m, 660],\n",
    "                    phi[m, 661], phi[m, 662], phi[m, 663], phi[m, 664], phi[m, 665], phi[m, 666], phi[m, 667], phi[m, 668], phi[m, 669], phi[m, 670],\n",
    "                    phi[m, 671], phi[m, 672], phi[m, 673], phi[m, 674], phi[m, 675], phi[m, 676], phi[m, 677], phi[m, 678], phi[m, 679], phi[m, 680],\n",
    "                    phi[m, 681], phi[m, 682], phi[m, 683], phi[m, 684], phi[m, 685], phi[m, 686], phi[m, 687], phi[m, 688], phi[m, 689], phi[m, 690],\n",
    "                    phi[m, 691], phi[m, 692], phi[m, 693], phi[m, 694], phi[m, 695], phi[m, 696], phi[m, 697], phi[m, 698], phi[m, 699], phi[m, 700],\n",
    "                    phi[m, 701], phi[m, 702], phi[m, 703], phi[m, 704], phi[m, 705], phi[m, 706], phi[m, 707], phi[m, 708], phi[m, 709], phi[m, 710],\n",
    "                    phi[m, 711], phi[m, 712], phi[m, 713], phi[m, 714], phi[m, 715], phi[m, 716], phi[m, 717], phi[m, 718], phi[m, 719], phi[m, 720],\n",
    "                    phi[m, 721], phi[m, 722], phi[m, 723], phi[m, 724], phi[m, 725], phi[m, 726], phi[m, 727], phi[m, 728], phi[m, 729], phi[m, 730],\n",
    "                    phi[m, 731], phi[m, 732], phi[m, 733], phi[m, 734], phi[m, 735], phi[m, 736], phi[m, 737], phi[m, 738], phi[m, 739], phi[m, 740],\n",
    "                    phi[m, 741], phi[m, 742], phi[m, 743], phi[m, 744], phi[m, 745], phi[m, 746], phi[m, 747], phi[m, 748], phi[m, 749], phi[m, 750],\n",
    "                    phi[m, 751], phi[m, 752], phi[m, 753], phi[m, 754], phi[m, 755], phi[m, 756], phi[m, 757], phi[m, 758], phi[m, 759], phi[m, 760],\n",
    "                    phi[m, 761], phi[m, 762], phi[m, 763], phi[m, 764], phi[m, 765], phi[m, 766], phi[m, 767], phi[m, 768], phi[m, 769], phi[m, 770],\n",
    "                    phi[m, 771], phi[m, 772], phi[m, 773], phi[m, 774], phi[m, 775], phi[m, 776], phi[m, 777], phi[m, 778], phi[m, 779], phi[m, 780],\n",
    "                    phi[m, 781], phi[m, 782], phi[m, 783], phi[m, 784], phi[m, 785], phi[m, 786], phi[m, 787], phi[m, 788], phi[m, 789], phi[m, 790],\n",
    "                    phi[m, 791], phi[m, 792], phi[m, 793], phi[m, 794], phi[m, 795], phi[m, 796], phi[m, 797], phi[m, 798], phi[m, 799], phi[m, 800],\n",
    "                    phi[m, 801], phi[m, 802], phi[m, 803], phi[m, 804], phi[m, 805], phi[m, 806], phi[m, 807], phi[m, 808], phi[m, 809], phi[m, 810],\n",
    "                    phi[m, 811], phi[m, 812], phi[m, 813], phi[m, 814], phi[m, 815], phi[m, 816], phi[m, 817], phi[m, 818], phi[m, 819], phi[m, 820],\n",
    "                    phi[m, 821], phi[m, 822], phi[m, 823], phi[m, 824], phi[m, 825], phi[m, 826], phi[m, 827], phi[m, 828], phi[m, 829], phi[m, 830],\n",
    "                    phi[m, 831], phi[m, 832], phi[m, 833], phi[m, 834], phi[m, 835], phi[m, 836], phi[m, 837], phi[m, 838], phi[m, 839], phi[m, 840],\n",
    "                    phi[m, 841], phi[m, 842], phi[m, 843], phi[m, 844], phi[m, 845], phi[m, 846], phi[m, 847], phi[m, 848], phi[m, 849], phi[m, 850],\n",
    "                    phi[m, 851], phi[m, 852], phi[m, 853], phi[m, 854], phi[m, 855], phi[m, 856], phi[m, 857], phi[m, 858], phi[m, 859], phi[m, 860],\n",
    "                    phi[m, 861], phi[m, 862], phi[m, 863], phi[m, 864], phi[m, 865], phi[m, 866], phi[m, 867], phi[m, 868], phi[m, 869], phi[m, 870],\n",
    "                    phi[m, 871], phi[m, 872], phi[m, 873], phi[m, 874], phi[m, 875], phi[m, 876], phi[m, 877], phi[m, 878], phi[m, 879], phi[m, 880],\n",
    "                    phi[m, 881], phi[m, 882], phi[m, 883], phi[m, 884], phi[m, 885], phi[m, 886], phi[m, 887], phi[m, 888], phi[m, 889], phi[m, 890],\n",
    "                    phi[m, 891], phi[m, 892], phi[m, 893], phi[m, 894], phi[m, 895], phi[m, 896], phi[m, 897], phi[m, 898], phi[m, 899], phi[m, 900]=#\n",
    "                )\n",
    "            for m = 1:length(βs)))\n",
    "    \n",
    "    if print_debug\n",
    "        println(\"JuMP model set up.\")\n",
    "    end\n",
    "    flush(stdout)\n",
    "    optimize!(model);\n",
    "    # println(\" ...solve status: $(status), mu = $(optimal_mu[t,i][j]), cont = $(optimal_cons[t,i][j])\")\n",
    "    return objective_value(model), [value(cons) value(mu)], termination_status(model)\n",
    "end\n",
    "\n",
    "\n",
    "function solve_timestep(stategrid, agents, t, Φ, Ψ, DICE_parameters, approx_parameters, algo_parameters, functions)\n",
    "    \n",
    "    # This function solves the problem for one entire timestep, agent by agent.\n",
    "    \n",
    "    @unpack initial_control, tries, verbosity = algo_parameters\n",
    "    @unpack a̲, a̅, n = approx_parameters\n",
    "    @unpack u, du, F, dF, Damage, dDamage, Abatecost, dAbatecost, K_next, S_next, eff, uneff = functions\n",
    "\n",
    "    #coefficients = Array{Array{Float64}}(undef, length(βs))\n",
    "\n",
    "    #global opt_cons, opt_mu, opt_value\n",
    "    opt_cons = Array{Array{Float64}}(undef, length(βs))\n",
    "    opt_mu = similar(opt_cons)\n",
    "    opt_value = similar(opt_cons)\n",
    "\n",
    "    for (i, agent) in enumerate(agents)                         # Solve the optimisation problem for each type\n",
    "        println(\"Agent $i\")\n",
    "\n",
    "        #opt_cons[i], opt_mu[i], opt_value[i] = solve_agent_problem(stategrid, i, t, Φ, Ψ, DICE_parameters, approx_parameters, algo_parameters, functions, [1e-8 1e-13])\n",
    "        opt_cons[i], opt_mu[i], opt_value[i] = solve_agent_problem(stategrid, i, t, Φ, Ψ, DICE_parameters, approx_parameters, algo_parameters, functions, [1e-8 1e-14])\n",
    "    \n",
    "   end\n",
    "\n",
    "    #return coefficients, V, opt_cons, opt_mu, opt_value\n",
    "    return [], opt_cons, opt_mu, opt_value\n",
    "    \n",
    "end\n",
    "\n",
    "\n",
    "\n",
    "function solve_agent_problem(stategrid, agent, t, Φ, Ψ, DICE_parameters, approx_parameters, algo_parameters, functions, tolerance)\n",
    "\n",
    "    # This function solves the agent's problem for a vector of states (typically the approximation nodes)\n",
    "    \n",
    "    @unpack initial_control, tries, verbosity = algo_parameters\n",
    "    @unpack a̲, a̅, n = approx_parameters\n",
    "    @unpack u, du, F, dF, Damage, dDamage, Abatecost, dAbatecost, K_next, S_next, eff, uneff = functions\n",
    "\n",
    "    opt_cons = zeros(size(stategrid, 1))\n",
    "    opt_mu = similar(opt_cons)\n",
    "    opt_value = similar(opt_cons)\n",
    "\n",
    "    for j in 1:size(stategrid, 1)                                  # Cycle through all the nodes\n",
    "        #println(\"State $j\")\n",
    "        current_Ke = stategrid[j, 1]\n",
    "        current_S = stategrid[j, 2]\n",
    "\n",
    "        initial_control = [uneff(current_Ke, t) .25]\n",
    "        if j > 1\n",
    "            if (j - 1) - n[1]*floor((j-1)/n[1]) == 0.0\n",
    "                initial_control = [opt_cons[j-n[1]] opt_mu[j - n[1]]]\n",
    "            else\n",
    "                initial_control = [opt_cons[j-1] opt_mu[j - 1]]\n",
    "            end\n",
    "        end\n",
    "\n",
    "        #println(\"Solving\")\n",
    "        model_results = solve_decision([current_Ke current_S], agent, t, Φ, DICE_parameters, approx_parameters, initial_control, functions, tolerance[1], verbosity)\n",
    "        attempt = 2;\n",
    "\n",
    "        while attempt <= tries\n",
    "\n",
    "            for i in 2:length(tolerance)\n",
    "                #print(\"...with tolerance $(tolerance[i])\")\n",
    "                model_results = solve_decision([current_Ke current_S], agent, t, Φ, DICE_parameters, approx_parameters, model_results[2], functions, tolerance[i], verbosity)\n",
    "            end\n",
    "            # @show typeof(model_results)\n",
    "            if ((model_results[3] == MOI.LOCALLY_SOLVED)  || (model_results[3] == MOI.NUMERICAL_ERROR))\n",
    "                break;\n",
    "            else\n",
    "                println(\"State $j (Ke = $(stategrid[j,1]), $(stategrid[j,2])) failed to solve on attempt $attempt.\")\n",
    "                println(model_results[3])\n",
    "                println(\"Non-optimal actions: c = $(model_results[2][1]), mu = $(model_results[2][2])\")\n",
    "                initial_control[2] = rand(1)[1]\n",
    "                initial_control[1] = rand(1)[1] * ((1 - δ) * uneff(current_Ke, t) / period_length + Y(A_seq[t], uneff(current_Ke, t), L_seq[t], initial_control[2], current_S, θ1[t]))\n",
    "                println(\"Trying again from c = $(initial_control[1]), mu = $(initial_control[2])\")\n",
    "                attempt = attempt + 1\n",
    "            end\n",
    "        end\n",
    "        # @show typeof(model_results)\n",
    "\n",
    "        # println(\"Solution\")\n",
    "        opt_value[j] = model_results[1]\n",
    "        opt_cons[j] = model_results[2][1]\n",
    "        opt_mu[j] = model_results[2][2]\n",
    "        status = model_results[3]\n",
    "\n",
    "        # println(\" ...solve status: $(status), mu = $(optimal_mu[t,i][j]), cont = $(optimal_cons[t,i][j])\")\n",
    "        @assert ((status == MOI.LOCALLY_SOLVED) || (status == MOI.NUMERICAL_ERROR))\n",
    "\n",
    "    end\n",
    "\n",
    "    return opt_cons, opt_mu, opt_value\n",
    "end\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# End of task 'initialise'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Task: load and/or solve the model.\n",
    "# To solve the model, the model must be first initialised (above).\n",
    "\n",
    "# This block used to load data for the currently specified run.\n",
    "# If loading a fully solved run, also construct the value functions by executing also the next block\n",
    "\n",
    "# If running from scratch, ignore this block and the next one; instead start solving frm period T.\n",
    "\n",
    "using JLD2, FileIO\n",
    "\n",
    "@load string(\".\\\\\", model_instance, \"\\\\\", solution_file, \".jld2\") stategrid Φ Values optimal_cons optimal_mu DICE_parameters approx_parameters functions #V\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reconstruct the functions V if required\n",
    "# Note: this needs to be done after loading the data.\n",
    "\n",
    "for (i, β) in enumerate(βs)\n",
    "    V[T, i] = (Ke, S) -> u.(((1 - δ) * uneff(Ke, T) + Y.(A_seq[T], uneff(Ke, T), L_seq[T], 0, S, θ1[T]))./L_seq[T], agents[i].η) * L_seq[T]\n",
    "end\n",
    "\n",
    "for t in T-1:-1:1                                              # Iterate backwards in time\n",
    "\n",
    "    for i in 1:length(agents)\n",
    "\n",
    "        V[t, i] = (Ke, S) -> (BasisMatrix(Basis(ChebParams(n[1], a̲[1], a̅[1]), ChebParams(n[2], a̲[2], a̅[2])), Expanded(), [Ke, S]).vals[1] * Φ[t,i])[1]\n",
    "        \n",
    "    end \n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now solve the outcome for period T, all types\n",
    "# This only needs to be done it period T has not been solved before\n",
    "#   (i.e. typically not if you are loading previously solved data)\n",
    "\n",
    "for (i, β) in enumerate(βs)\n",
    "    V[T, i] = (Ke, S) -> u.(((1 - δ) * uneff(Ke, T) + Y.(A_seq[T], uneff(Ke, T), L_seq[T], 0, S, θ1[T]))./L_seq[T], agents[i].η) * L_seq[T]\n",
    "    optimal_cons[T,i] = (1 - δ) * uneff(stategrid[:, 1], T) + Y.(A_seq[T], uneff(stategrid[:, 1], T), L_seq[T], 0, stategrid[:, 2], θ1[T])\n",
    "    optimal_mu[T,i] = zeros(size(stategrid, 1))\n",
    "    Values[T,i] = V[T, i].(stategrid[:,1], stategrid[:,2])\n",
    "    \n",
    "    Φ[T,i] = Ψ \\ Values[T,i]\n",
    "\n",
    "end\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now solve for periods t < T.\n",
    "# Again, this only needs to be done once, after which the data can be saved.\n",
    "\n",
    "# Model solution may take a while and I have found it safer to run a few timesteps at a time, then use the debug functions\n",
    "#   (see below) to ensure the value functions remain well-behaved.\n",
    "\n",
    "for t in T-1:-1:1                                             # Iterate backwards in time\n",
    "    println(\"Timestep $t\")\n",
    "\n",
    "    timestep_solution = solve_timestep(stategrid, agents, t, Φ, Ψ, DICE_parameters, approx_parameters, algo_parameters, functions)\n",
    "    for i in 1:length(agents)\n",
    "        #V[t,i] = timestep_solution[2][i]\n",
    "        optimal_cons[t,i] = timestep_solution[2][i]\n",
    "        optimal_mu[t,i] = timestep_solution[3][i]\n",
    "        Values[t,i] = timestep_solution[4][i]\n",
    "        \n",
    "\n",
    "        Φ[t,i] = Ψ \\ Values[t,i]\n",
    "\n",
    "        V[t, i] = (Ke, S) -> (BasisMatrix(Basis(ChebParams(n[1], a̲[1], a̅[1]), ChebParams(n[2], a̲[2], a̅[2])), Expanded(), [Ke, S]).vals[1] * Φ[t,i])[1]\n",
    "\n",
    "    end \n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Code for saving the solution computed thus far.\n",
    "\n",
    "using JLD2, FileIO\n",
    "\n",
    "@save string(\".\\\\\", model_instance, \"\\\\\", solution_file, \".jld2\") stategrid Φ Values optimal_cons optimal_mu DICE_parameters approx_parameters functions V"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# End of task 'solve the model'."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Task: debug or construct approximation errors.\n",
    "\n",
    "# This block for debugging purposes only.  Computes the maximand at a given (state, time, agent) for some specified\n",
    "#   range of values of the choice variables (consumption, abatement rate)\n",
    "\n",
    "# Note: use this block only if you need to diagnose problems in the optimisation problem for a given time, type, and state.\n",
    "# The next two block is more useful for detecting problems in the value function iteration process.\n",
    "\n",
    "using Plots\n",
    "using Statistics\n",
    "\n",
    "test_t = 45\n",
    "test_agent = 20\n",
    "\n",
    "test_node = 192\n",
    "\n",
    "# Different options for specifying the state and the choice variables around which the range is constructed\n",
    "\n",
    "#c_test_K = test_grid[test_node][1]\n",
    "#c_test_S = test_grid[test_node][2]\n",
    "#c_test_cons = test_cons[test_node]\n",
    "#c_test_mu = test_mu[test_node]\n",
    "\n",
    "c_test_K = stategrid[test_node, 1]\n",
    "c_test_S = stategrid[test_node, 2]\n",
    "c_test_cons = 7376.74858\n",
    "c_test_mu = .999999\n",
    "\n",
    "#c_test_cons = optimal_cons[test_t, test_agent][test_node]\n",
    "#c_test_mu = optimal_mu[test_t, test_agent][test_node]\n",
    "#c_test_cons = 879\n",
    "#c_test_mu = .7\n",
    "\n",
    "println(\"Testing optimality of decision: period = $test_t, agent = $test_agent, K = $(c_test_K), S = $(c_test_S)\")\n",
    "\n",
    "# The desired ranges and resolution can be modified below\n",
    "c_test_mu_range = range(0.99*c_test_mu, stop = min(1.0, 1.01*c_test_mu), length=100)\n",
    "c_test_cons_range = range(0.99*c_test_cons, stop = 1.01*c_test_cons, length=100)\n",
    "\n",
    "test_val = (mu, cons) -> u.(cons./L_seq[test_t], ηs[test_agent]) * L_seq[test_t] .+ βs[test_agent] .* sum(weights[test_agent, m] .* V[test_t + 1, m].(eff((1 - δ) .* uneff(c_test_K, test_t) .+ period_length*(Y.(A_seq[test_t], uneff(c_test_K, test_1), L_seq[test_t], mu, c_test_S, θ1[test_t]) .- cons), test_t + 1), c_test_S .+ period_length*(σ[test_t] .* (1 .- mu) .* F(A_seq[test_t], uneff(c_test_K, test_t), L_seq[test_t]) + E_land[test_t])/1000) for m = 1:length(βs));\n",
    "\n",
    "test_val2 = (mu, cons) -> u.(cons./L_seq[test_t], ηs[test_agent])*L_seq[test_t] .+ βs[test_agent] .* sum(weights[test_agent, m] .* (BasisMatrix(basis, Expanded(), [eff((1 - δ) .* uneff(c_test_K, test_t) .+ period_length*(Y.(A_seq[test_t], uneff(c_test_K, test_t), L_seq[test_t], mu, c_test_S, θ1[test_t]) .- cons), test_t + 1),\n",
    "                c_test_S .+ period_length*(σ[test_t] .* (1 .- mu) .* F(A_seq[test_t], uneff(c_test_K, test_t), L_seq[test_t]) + E_land_seq[test_t])/1000]).vals[1] * Φ[test_t+1,m]) for m = 1:length(βs))\n",
    "\n",
    "test_val3 = (mu, cons) -> βs[test_agent] .* sum(weights[test_agent, m] .* (BasisMatrix(basis, Expanded(), [eff((1 - δ) .* uneff(c_test_K, test_t) .+ period_length*(Y.(A_seq[test_t], uneff(c_test_K, test_t), L_seq[test_t], mu, c_test_S, θ1[test_t]) .- cons), test_t + 1),\n",
    "                c_test_S .+ period_length*(σ[test_t] .* (1 .- mu) .* F(A_seq[test_t], uneff(c_test_K, test_t), L_seq[test_t]) + E_land_seq[test_t])/1000]).vals[1] * Φ[test_t+1,m]) for m = 1:length(βs))\n",
    "\n",
    "test_val4 = (mu, cons) -> u.(cons./L_seq[test_t], ηs[test_agent])*L_seq[test_t]\n",
    "\n",
    "test_Ke_next = (mu, cons) -> eff((1 - δ) .* uneff(c_test_K, test_t) .+ period_length*(Y.(A_seq[test_t], uneff(c_test_K, test_t), L_seq[test_t], mu, c_test_S, θ1[test_t]) .- cons), test_t + 1)\n",
    "test_Se_next = (mu, cons) -> c_test_S .+ period_length*(σ[test_t] .* (1 .- mu) .* F(A_seq[test_t], uneff(c_test_K, test_t), L_seq[test_t]) + E_land_seq[test_t]) / 1000\n",
    "\n",
    "plt1 = plot(c_test_mu_range, title=\"Plot 1\", [a[1] for a in test_val2.(c_test_mu_range, c_test_cons)]);\n",
    "vline!([c_test_mu], color=[:red], width=[1.], alpha=[0.9]);\n",
    "display(plt1)\n",
    "\n",
    "plt2 = plot(c_test_cons_range, title=\"Plot 2\", [a[1] for a in test_val2.(c_test_mu, c_test_cons_range)], reuse=false)\n",
    "vline!([c_test_cons], color=[:red], width=[1.], alpha=[0.9])\n",
    "display(plt2)\n",
    "\n",
    "plt5 = plot(c_test_cons_range, title=\"Plot 5\", [a[1] for a in test_Ke_next.(c_test_mu, c_test_cons_range)], reuse=false)\n",
    "display(plt5)\n",
    "\n",
    "plt7 = plot(c_test_cons_range, title=\"Plot 7\", [a[1] for a in V[test_t + 1, test_agent].(eff.(test_Ke_next.(c_test_mu, c_test_cons_range), test_t + 1), test_Se_next.(c_test_mu, c_test_cons_range))], reuse=false)\n",
    "display(plt7)\n",
    "\n",
    "plt6 = plot(c_test_mu_range, title=\"Plot 6\", [a[1] for a in test_Se_next.(c_test_mu_range, c_test_cons)], reuse=false)\n",
    "display(plt6)\n",
    "\n",
    "plt3 = plot(c_test_cons_range, title=\"Plot 3\", [a[1] for a in test_val3.(c_test_mu, c_test_cons_range)], reuse=false)\n",
    "vline!([c_test_cons], color=[:blue], width=[1.], alpha=[0.9])\n",
    "\n",
    "plt4 = plot(c_test_cons_range, title=\"Plot 4\", [a[1] for a in test_val4.(c_test_mu, c_test_cons_range)], reuse=false)\n",
    "vline!([c_test_cons], color=[:red], width=[1.], alpha=[0.9])\n",
    "\n",
    "display(plt1)\n",
    "display(plt2)\n",
    "display(plt3)\n",
    "display(plt4)\n",
    "\n",
    "println(\"The implied optimum yields next period K = $(test_Ke_next(c_test_mu, c_test_cons)), S = $(test_Se_next(c_test_mu, c_test_cons)).\");\n",
    "println(\"Compare these with $a̲[1] < K < $a̅[1], $a̲[2] < S < $a̅[2].\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Debugging and robustness check purposes only.\n",
    "# Solves the value function for a given time, agent at some state grid (typically NOT the approximation grid)\n",
    "#   and computes various relative errors and derivatives of the value function, and the optimal choice variables.\n",
    "# Next block is then used to plot these. Problems are typically easily identified by visual inspection.\n",
    "\n",
    "# But this is an important check on the numerical stability. Instability is pretty apparent is the figures\n",
    "#   produced and usually appears at high levels of K.  Solution is to increase the range of approximation\n",
    "#   (to higher values of effective capital) and start the solution over from T down.\n",
    "\n",
    "using Plots\n",
    "\n",
    "test_time = 1\n",
    "test_agent = 1\n",
    "\n",
    "test_Kegrid = range(1.05*Kegrid[1], .95 * Kegrid[end], length=37)\n",
    "test_Sgrid = range(1.05*Sgrid[1], .95 * Sgrid[end], length=37)\n",
    "test_grid = vec(collect(Base.product(test_Kegrid, test_Sgrid)))\n",
    "test_value = zeros(size(test_grid, 1))\n",
    "test_value_exact = similar(test_value)\n",
    "test_error = similar(test_value)\n",
    "test_cons = similar(test_value)\n",
    "test_mu = similar(test_value)\n",
    "test_error_K = similar(test_value)\n",
    "test_value_K = similar(test_value)\n",
    "test_value_KK = similar(test_value)\n",
    "test_value_SS = similar(test_value)\n",
    "test_value_SK = similar(test_value)\n",
    "test_value_S = similar(test_value)\n",
    "test_value_K_exact = similar(test_value)\n",
    "test_du = similar(test_value)\n",
    "\n",
    "\n",
    "attempt = 1;\n",
    "tries = 12;\n",
    "initial_control = [600.0, .5]\n",
    "\n",
    "for j in 1:size(test_grid,1)\n",
    "#for j in 1002:1002\n",
    "#for j in 48:size(test_grid,1)\n",
    "    println(\"State $j: K = $(test_grid[j][1]), S = $(test_grid[j][2])\")\n",
    "    attempt = 1\n",
    "    if test_time < T\n",
    "        global test_K_index = 1\n",
    "        global test_S_index = 1\n",
    "        \n",
    "        for (i, h) in enumerate(collect(Kegrid))\n",
    "            if h > test_grid[j][1]\n",
    "                #print(i)\n",
    "                test_K_index = i;\n",
    "                break;\n",
    "            else\n",
    "                #print(i)\n",
    "            end\n",
    "        end\n",
    "        for (i, h) in enumerate(collect(Sgrid))\n",
    "            if h > test_grid[j][2]\n",
    "                #print(i)\n",
    "                test_S_index = i;\n",
    "                break;\n",
    "            else\n",
    "                #print(i)\n",
    "            end\n",
    "        end\n",
    "        initial_control = (1 - (test_grid[j][1] - Kegrid[test_K_index-1])/(Kegrid[test_K_index] - Kegrid[test_K_index-1])) * (1 - (test_grid[j][2] - Sgrid[test_S_index-1])/(Sgrid[test_S_index] - Sgrid[test_S_index-1])) *\n",
    "                                           [optimal_cons[test_time, test_agent][(test_S_index - 2) * n[1] + test_K_index] optimal_mu[test_time, test_agent][(test_S_index - 2) * n[1] + test_K_index]] +\n",
    "                                    (test_grid[j][1] - Kegrid[test_K_index-1])/(Kegrid[test_K_index] - Kegrid[test_K_index-1]) * (1 - (test_grid[j][2] - Sgrid[test_S_index-1])/(Sgrid[test_S_index] - Sgrid[test_S_index-1])) *\n",
    "                                           [optimal_cons[test_time, test_agent][(test_S_index - 2) * n[1] + test_K_index + 1] optimal_mu[test_time, test_agent][(test_S_index - 2) * n[1] + test_K_index + 1]] +\n",
    "                                    (1 - (test_grid[j][1] - Kegrid[test_K_index-1])/(Kegrid[test_K_index] - Kegrid[test_K_index-1])) * (test_grid[j][2] - Sgrid[test_S_index-1])/(Sgrid[test_S_index] - Sgrid[test_S_index-1]) *\n",
    "                                           [optimal_cons[test_time, test_agent][(test_S_index - 1) * n[1] + test_K_index] optimal_mu[test_time, test_agent][(test_S_index - 1) * n[1] + test_K_index]] +\n",
    "                                    (test_grid[j][1] - Kegrid[test_K_index-1])/(Kegrid[test_K_index] - Kegrid[test_K_index-1]) * (test_grid[j][2] - Sgrid[test_S_index-1])/(Sgrid[test_S_index] - Sgrid[test_S_index-1]) *\n",
    "                                           [optimal_cons[test_time, test_agent][(test_S_index - 1) * n[1] + test_K_index + 1] optimal_mu[test_time, test_agent][(test_S_index - 1) * n[1] + test_K_index + 1]]\n",
    "        #= println(\"Initial guess: ($(initial_control[1]), $(initial_control[2])), weighted of ($(optimal_cons[test_time, test_agent][(test_S_index - 2) * n[1] + test_K_index]), $(optimal_mu[test_time, test_agent][(test_S_index - 2) * n[1] + test_K_index])),\n",
    "                    ($(optimal_cons[test_time, test_agent][(test_S_index - 2) * n[1] + test_K_index + 1]), $(optimal_mu[test_time, test_agent][(test_S_index - 2) * n[1] + test_K_index + 1])),\n",
    "                    ($(optimal_cons[test_time, test_agent][(test_S_index - 1) * n[1] + test_K_index]), $(optimal_mu[test_time, test_agent][(test_S_index - 1) * n[1] + test_K_index])),\n",
    "                    ($(optimal_cons[test_time, test_agent][(test_S_index - 1) * n[1] + test_K_index + 1]), $(optimal_mu[test_time, test_agent][(test_S_index - 1) * n[1] + test_K_index + 1]))\");\n",
    "        =#\n",
    "        \n",
    "        while(attempt <= tries)\n",
    "            # print(\"1\")\n",
    "            global test_sol = solve_decision(collect(test_grid[j, :][1]), test_agent, test_time, Φ, DICE_parameters, approx_parameters, initial_control, functions, 1e-8, 0)\n",
    "            global test_sol = solve_decision(collect(test_grid[j, :][1]), test_agent, test_time, Φ, DICE_parameters, approx_parameters, test_sol[2], functions, 1e-14, 0)\n",
    "            if ((test_sol[3] == MOI.LOCALLY_SOLVED) || (test_sol[3] == MOI.NUMERICAL_ERROR))\n",
    "                println(\"Status: $(test_sol[3])\")\n",
    "                break;\n",
    "            else \n",
    "                println(\"State $j (K = $(test_grid[j][1]), $(test_grid[j][2])) failed to solve on attempt $attempt.\")\n",
    "                println(test_sol[3])\n",
    "                initial_control[2] = rand(1)[1]\n",
    "                initial_control[1] = rand(1)[1] * eff(((1 - δ) * uneff(test_grid[j][1], test_time) / period_length + Y(A_seq[test_time], uneff(test_grid[j][1], test_time), L_seq[test_time], initial_control[2], test_grid[j][2], θ1[test_time])), test_time)\n",
    "                println(\"Trying again from c = $initial_control[1], mu = $initial_control[2]\")\n",
    "                attempt = attempt + 1\n",
    "            end\n",
    "        end\n",
    "        test_value[j] = test_sol[1]\n",
    "        test_error[j] = test_value[j] - (BasisMatrix(Basis(ChebParams(n[1], a̲[1], a̅[1]), ChebParams(n[2], a̲[2], a̅[2])), \n",
    "            Expanded(), collect(test_grid[j,:][1])).vals[1] * collect(Φ[test_time, test_agent]))[1]\n",
    "        test_cons[j] = test_sol[2][1]\n",
    "        test_mu[j] = test_sol[2][2]        \n",
    "        test_value_K[j] = (BasisMatrix(Basis(ChebParams(n[1], a̲[1], a̅[1]), ChebParams(n[2], a̲[2], a̅[2])), \n",
    "            Expanded(), collect(test_grid[j,:][1]), [1 0]).vals[1] * collect(Φ[test_time, test_agent]))[1]\n",
    "        test_value_KK[j] = (BasisMatrix(Basis(ChebParams(n[1], a̲[1], a̅[1]), ChebParams(n[2], a̲[2], a̅[2])), \n",
    "            Expanded(), collect(test_grid[j,:][1]), [2 0]).vals[1] * collect(Φ[test_time, test_agent]))[1]\n",
    "        test_value_SS[j] = (BasisMatrix(Basis(ChebParams(n[1], a̲[1], a̅[1]), ChebParams(n[2], a̲[2], a̅[2])), \n",
    "            Expanded(), collect(test_grid[j,:][1]), [0 2]).vals[1] * collect(Φ[test_time, test_agent]))[1]\n",
    "        test_value_SK[j] = (BasisMatrix(Basis(ChebParams(n[1], a̲[1], a̅[1]), ChebParams(n[2], a̲[2], a̅[2])), \n",
    "            Expanded(), collect(test_grid[j,:][1]), [1 1]).vals[1] * collect(Φ[test_time, test_agent]))[1]\n",
    "        test_value_S[j] = (BasisMatrix(Basis(ChebParams(n[1], a̲[1], a̅[1]), ChebParams(n[2], a̲[2], a̅[2])), \n",
    "            Expanded(), collect(test_grid[j,:][1]), [0 1]).vals[1] * collect(Φ[test_time, test_agent]))[1]\n",
    "        test_value_K_exact[j] = agents[test_agent].β * sum(weights[test_agent, m] *\n",
    "                (BasisMatrix(Basis(ChebParams(n[1], a̲[1], a̅[1]), ChebParams(n[2], a̲[2], a̅[2])), \n",
    "                    Expanded(), [eff((1 - δ) * uneff(test_grid[j][1], test_time) + period_length*(Y.(A_seq[test_time], uneff(test_grid[j][1], test_time), L_seq[test_time], test_mu[j], test_grid[j][2], θ1[test_time]) - test_cons[j]), test_time + 1), test_grid[j][2] + period_length * (σ[test_time] * (1 - test_mu[j]) * F(A_seq[test_time], uneff(test_grid[j][1], test_time), L_seq[test_time]) + E_land_seq[test_time])/ 1000], [1 0]).vals[1] * collect(Φ[test_time + 1, m]))[1] *\n",
    "                eff(1, test_time + 1) * (1 - δ + period_length * dY.(A_seq[test_time], uneff(test_grid[j][1], test_time), L_seq[test_time], test_mu[j], test_grid[j][2], θ1[test_time])) * uneff(1, test_time) +\n",
    "                (BasisMatrix(Basis(ChebParams(n[1], a̲[1], a̅[1]), ChebParams(n[2], a̲[2], a̅[2])), \n",
    "                    Expanded(), [eff((1 - δ) * uneff(test_grid[j][1], test_time) + period_length * (Y.(A_seq[test_time], uneff(test_grid[j][1], test_time), L_seq[test_time], test_mu[j], test_grid[j][2], θ1[test_time]) - test_cons[j]), test_time + 1), test_grid[j][2] + period_length * (σ[test_time] * (1 - test_mu[j]) * F(A_seq[test_time], uneff(test_grid[j][1], test_time), L_seq[test_time]) + E_land_seq[test_time]) / 1000], [0 1]).vals[1] * collect(Φ[test_time + 1, m]))[1] *\n",
    "                (period_length * σ[test_time] * (1 - test_mu[j]) * dF(A_seq[test_time], uneff(test_grid[j][1], test_time), L_seq[test_time]) / 1000)\n",
    "            for m in 1:length(βs))\n",
    "        test_error_K[j] = test_value_K_exact[j] - test_value_K[j]\n",
    "    else\n",
    "        test_value_exact[j] = u(((1 - δ) * uneff(test_grid[j][1], T) + Y.(A_seq[T], uneff(test_grid[j][1], T), L_seq[T], 0, test_grid[j][2], θ1[T]))./L_seq[T], agents[test_agent].η) * L_seq[T]\n",
    "        test_value[j] = (BasisMatrix(Basis(ChebParams(n[1], a̲[1], a̅[1]), ChebParams(n[2], a̲[2], a̅[2])), \n",
    "            Expanded(), collect(test_grid[j,:][1])).vals[1] * collect(Φ[test_time, test_agent]))[1]\n",
    "        test_error[j] = test_value_exact[j] - test_value[j]\n",
    "        test_cons[j] = ((1 - δ) * uneff(test_grid[j][1], T) + Y.(A_seq[T], uneff(test_grid[j][1], T), L_seq[T], 0, test_grid[j][2], θ1[T]));\n",
    "        test_mu[j] = 0;\n",
    "        test_value_K[j] = (BasisMatrix(Basis(ChebParams(n[1], a̲[1], a̅[1]), ChebParams(n[2], a̲[2], a̅[2])), \n",
    "            Expanded(), collect(test_grid[j,:][1]), [1 0]).vals[1] * collect(Φ[test_time, test_agent]))[1]\n",
    "        test_value_S[j] = (BasisMatrix(Basis(ChebParams(n[1], a̲[1], a̅[1]), ChebParams(n[2], a̲[2], a̅[2])), \n",
    "            Expanded(), collect(test_grid[j,:][1]), [0 1]).vals[1] * collect(Φ[test_time, test_agent]))[1]\n",
    "        test_value_KK[j] = (BasisMatrix(Basis(ChebParams(n[1], a̲[1], a̅[1]), ChebParams(n[2], a̲[2], a̅[2])), \n",
    "            Expanded(), collect(test_grid[j,:][1]), [2 0]).vals[1] * collect(Φ[test_time, test_agent]))[1]\n",
    "        test_value_K_exact[j] = uneff(du(((1 - δ) * uneff(test_grid[j][1], T) + Y.(A_seq[T], uneff(test_grid[j][1], T), L_seq[T], 0, test_grid[j][2], θ1[T]))./L_seq[T], agents[test_agent].η) * ((1 - δ) + dY.(A_seq[T], uneff(test_grid[j][1], T), L_seq[T], 0, test_grid[j][2], θ1[T])), T)\n",
    "        test_error_K[j] = test_value_K_exact[j] - test_value_K[j]\n",
    "        test_du[j] = du(((1 - δ) * uneff(test_grid[j][1], T) + Y.(A_seq[T], uneff(test_grid[j][1], T), L_seq[T], 0, test_grid[j][2], θ1[T]))./L_seq[T], agents[test_agent].η)\n",
    "    end\n",
    "\n",
    "    \n",
    "end\n",
    "\n",
    "maximal_error = findmax(abs.(test_error./test_value))\n",
    "maximal_error_V_K = findmax(abs.(test_error_K./test_value_K))\n",
    "\n",
    "println(\"The maximal relative error for agent = $(test_agent) at time t = $(test_time) is $(maximal_error[1]) at K = $(test_grid[maximal_error[2]][1]), S = $(test_grid[maximal_error[2]][2]); for V_K = $(maximal_error_V_K[1])\")\n",
    "\n",
    "Vplot = (K = zeros(length(test_grid)), S = zeros(length(test_grid)), V = zeros(length(test_grid)), error = zeros(length(test_grid)), cons = zeros(length(test_grid)), mu = zeros(length(test_grid)), V_K = zeros(length(test_grid)), V_S = zeros(length(test_grid)), V_KK = zeros(length(test_grid)), V_SK = zeros(length(test_grid)), V_SS = zeros(length(test_grid)), error_V_K = zeros(length(test_grid)), error_V_K_abs = zeros(length(test_grid)), du = zeros(length(test_grid)))\n",
    "\n",
    "for i = 1:length(test_grid)\n",
    "    Vplot.K[i] = test_grid[i][1]\n",
    "    Vplot.S[i] = test_grid[i][2]\n",
    "    Vplot.V[i] = test_value[i]\n",
    "    Vplot.error[i] = test_error[i]./test_value[i]\n",
    "    Vplot.cons[i] = test_cons[i]\n",
    "    Vplot.mu[i] = test_mu[i]\n",
    "    Vplot.V_K[i] = test_value_K[i]\n",
    "    Vplot.V_S[i] = test_value_S[i]\n",
    "    Vplot.V_KK[i] = test_value_KK[i]\n",
    "    Vplot.V_SS[i] = test_value_SS[i]\n",
    "    Vplot.V_SK[i] = test_value_SK[i]\n",
    "    Vplot.error_V_K[i] = test_error_K[i] ./ test_value_K[i]\n",
    "    Vplot.error_V_K_abs[i] = test_error_K[i]\n",
    "    Vplot.du[i] = test_du[i]\n",
    "end\n",
    "    \n",
    "pyplot()\n",
    "plot(Vplot.K, Vplot.S, Vplot.error, st=:surface, camera=(20,40))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plots the value function, some derivatives, and the optimal control variables as function of state.\n",
    "# These were computed in the previous block.\n",
    "# Numerical instability is typically obvious by visual inspection, although make sure that the z-scale is\n",
    "#   of meaningful magnitudes.\n",
    "\n",
    "plt1 = plot(Vplot.K, Vplot.S, Vplot.V, st=:surface, camera=(-45,30));\n",
    "display(plt1)\n",
    "\n",
    "plt2 = plot(Vplot.K, Vplot.S, Vplot.V_K, st=:surface, camera=(10,30))\n",
    "display(plt2)\n",
    "\n",
    "plt3 = plot(Vplot.K, Vplot.S, test_value_KK, st=:surface, camera=(20,20))\n",
    "display(plt3)\n",
    "\n",
    "plt4 = plot(Vplot.K, Vplot.S, Vplot.error_V_K, st=:surface, camera=(-30,10))\n",
    "display(plt4)\n",
    "\n",
    "plt5 = plot(Vplot.K, Vplot.S, Vplot.cons, st=:surface, camera=(-35,15))\n",
    "display(plt5)\n",
    "\n",
    "plt6 = plot(Vplot.K, Vplot.S, Vplot.mu, st=:surface, camera=(20,20));\n",
    "display(plt6)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# End of task: debuggin functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Task: initialise tools for analysing results.\n",
    "# These functions are used a lot for various purposes below.\n",
    "\n",
    "# This function is required to compute trajectories of the system. Used to plot figures, run Monte Carlo, etc.\n",
    "\n",
    "function get_trajectories(initial_state, Φ, horizon, type_seq, DICE_parameters, approx_parameters, functions)\n",
    "    \n",
    "    @unpack A, L_seq, ηs, βs, weights, θ1, δ, σ, T, Ae, E_land = DICE_parameters\n",
    "    @unpack a̲, a̅, n, period_length = approx_parameters\n",
    "    @unpack u, du, F, dF, Damage, dDamage, Abatecost, dAbatecost, K_next, S_next, eff, uneff = functions\n",
    "    \n",
    "    K_next(K, Y, C) = (1 - δ) * K + period_length*(Y - C)\n",
    "    S_next(S, F, μ, σ, E_land) = S + period_length*(σ * (1 - μ) * F + E_land) / 1000\n",
    "\n",
    "    initial_K = initial_state[1]\n",
    "    initial_S = initial_state[2]\n",
    "\n",
    "    # types = [1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1]\n",
    "    \n",
    "    paths = (K_t = zeros(T), S_t = zeros(T), Y_t = zeros(T), F_t = zeros(T), c_t = zeros(T), μ_t = zeros(T), u_t = zeros(T), E_t = zeros(T), p_t = zeros(T), SCC_t = zeros(T));\n",
    "    #paths.K_t = ones(T)\n",
    "    #paths.S_t = similar(K_t)\n",
    "    #paths.Y_t = similar(K_t)\n",
    "    #F_t = similar(K_t)\n",
    "    #c_t = similar(K_t)\n",
    "    #μ_t = similar(K_t)\n",
    "    #u_t = similar(K_t)\n",
    "    #E_t = similar(K_t)\n",
    "    #p_t = similar(K_t)\n",
    "\n",
    "    paths.K_t[1] = initial_K\n",
    "    paths.S_t[1] = initial_S\n",
    "    soln = solve_decision([paths.K_t[1] paths.S_t[1]], type_seq[1], 1, Φ, DICE_parameters, approx_parameters, [uneff(paths.K_t[1]/2, 1) .3], functions, 1e-6, 0)\n",
    "    soln = solve_decision([paths.K_t[1] paths.S_t[1]], type_seq[1], 1, Φ, DICE_parameters, approx_parameters, soln[2], functions, 1e-10, 0)\n",
    "    paths.c_t[1] = soln[2][1]\n",
    "    paths.μ_t[1] = soln[2][2]\n",
    "    paths.Y_t[1] = Y(A[1], uneff(paths.K_t[1], 1), L_seq[1], paths.μ_t[1], paths.S_t[1], θ1[1])\n",
    "    paths.F_t[1] = F(A[1], uneff(paths.K_t[1], 1), L_seq[1])\n",
    "    paths.E_t[1] = σ[1] * (1 - paths.μ_t[1]) * paths.F_t[1]\n",
    "    paths.u_t[1] = u(paths.c_t[1]/L_seq[1], ηs[type_seq[1]])\n",
    "\n",
    "    for t = 1:min(horizon, T-1)\n",
    "        soln = solve_decision([paths.K_t[t] paths.S_t[t]], type_seq[t], t, Φ, DICE_parameters, approx_parameters, [uneff(paths.K_t[t]/2, t) .3], functions, 1e-6, 0)\n",
    "        soln = solve_decision([paths.K_t[t] paths.S_t[t]], type_seq[t], t, Φ, DICE_parameters, approx_parameters, soln[2], functions, 1e-10, 0)\n",
    "        paths.c_t[t] = soln[2][1]\n",
    "        paths.μ_t[t] = soln[2][2]\n",
    "        paths.Y_t[t] = Y(A[t], uneff(paths.K_t[t], t), L_seq[t], paths.μ_t[t], paths.S_t[t], θ1[t])\n",
    "        paths.F_t[t] = F(A[t], uneff(paths.K_t[t], t), L_seq[t])\n",
    "        paths.u_t[t] = u(paths.c_t[t]/L_seq[t], ηs[type_seq[t]])\n",
    "        paths.E_t[t] = σ[t] * (1 - paths.μ_t[t]) * paths.F_t[t]    # E is just the industrial emissions!\n",
    "        paths.K_t[t+1] = eff(K_next(uneff(paths.K_t[t], t), paths.Y_t[t], paths.c_t[t]), t + 1)\n",
    "        paths.S_t[t+1] = S_next(paths.S_t[t], paths.F_t[t], paths.μ_t[t], σ[t], E_land[t])\n",
    "        paths.SCC_t[t] = -sum(weights[type_seq[t], m] * (BasisMatrix(Basis(ChebParams(n[1], a̲[1], a̅[1]), ChebParams(n[2], a̲[2], a̅[2])), \n",
    "                       Expanded(), [paths.K_t[t + 1] paths.S_t[t + 1]], [0 1]).vals[1] * collect(Φ[t+1, m]))[1] for m in 1:length(βs)) /\n",
    "                   sum(weights[type_seq[t], m] * (BasisMatrix(Basis(ChebParams(n[1], a̲[1], a̅[1]), ChebParams(n[2], a̲[2], a̅[2])), \n",
    "                       Expanded(), [paths.K_t[t + 1] paths.S_t[t + 1]], [1 0]).vals[1] * collect(Φ[t+1, m]))[1] for m in 1:length(βs)) *\n",
    "                   uneff(1, t + 1)\n",
    "        paths.p_t[t] = dAbatecost(paths.μ_t[t], θ1[t]) / (σ[t] / 1000); # Should the 1000 not be multiplied in the denominator? Adjusts for tC, from $tn/GtC...\n",
    "    end\n",
    "\n",
    "    if (horizon >= T)\n",
    "        paths.p_t[T] = 0.0\n",
    "        paths.μ_t[T] = 0\n",
    "        paths.Y_t[T] = Y(A[T], uneff(paths.K_t[T], T), L_seq[T], paths.μ_t[T], paths.S_t[T], θ1[T])\n",
    "        paths.c_t[T] = paths.Y_t[T] + (1 - δ) * uneff(paths.K_t[T], T)\n",
    "        paths.F_t[T] = F(A[T], uneff(paths.K_t[T], T), L_seq[T])\n",
    "        paths.u_t[T]  = u(paths.c_t[T]/L_seq[T], ηs[type_seq[T]])\n",
    "        paths.E_t[T] = σ[T] * (1 - paths.μ_t[T]) * paths.F_t[T]\n",
    "    end\n",
    "    \n",
    "    return paths\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to get a set of matrices used in plotting trajectories of different types\n",
    "\n",
    "function get_plotdata(paths, T, n)\n",
    "\n",
    "    plotdata = (c = zeros(n, T), K = zeros(n, T), S = zeros(n, T), p = zeros(n, T), μ = zeros(n, T), SCC = zeros(n, T))\n",
    "    \n",
    "    for j = 1:n\n",
    "        plotdata.c[j, :] = paths[j].c_t\n",
    "        plotdata.K[j, :] = paths[j].K_t\n",
    "        plotdata.S[j, :] = paths[j].S_t\n",
    "        plotdata.p[j, :] = paths[j].p_t    \n",
    "        plotdata.SCC[j, :] = paths[j].SCC_t    \n",
    "        plotdata.μ[j, :] = paths[j].μ_t    \n",
    "    end\n",
    "    \n",
    "    return plotdata\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# End of task: initialise functions for analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Task: Compute the DICE benchmark\n",
    "# This task only used to produce the 'DICE' plot looking at nondogmatic SCC values computed at various times along the\n",
    "#   standard DICE trajectory.\n",
    "# Note: used in conjunction with the next task. This requires:\n",
    "#   (1) Solving the non-dogmatic model fully and saving the data.\n",
    "#   (2) Reinitialising for the DICE/Stern benchmark parameterisation\n",
    "#   (3) Computing the DICE path (which is deterministic) using the next block. This produces a variable called 'paths'.\n",
    "#   (4) Using this block, renaming the 'paths' variable and saving it.\n",
    "#   SEE BELOW\n",
    "\n",
    "# This function useful for saving the benchmark data\n",
    "\n",
    "paths_benchmark = paths\n",
    "@save string(\".\\\\\", model_instance, \"\\\\paths_\", model_instance, \"Benchmark.jld2\") paths_benchmark\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use this to load paths_benchmark if needed.\n",
    "\n",
    "@load string(\".\\\\Y250T10N2Benchmark\\\\paths_Y250T10N2Benchmark.jld2\") paths_benchmark"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialise the structure to save the SCCs computed starting at various points along the DICE benchmark trajectory\n",
    "\n",
    "# First parameter: how many non-dogmatic agents\n",
    "# Second parameter: how many different degrees of nondogmatism (Figure 4 has three of these)\n",
    "# Third parameter: how may initial times (Figure 4 has five of these)\n",
    "DICE_SCC_pairs = zeros(agent_type_n, 3, 5)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This computes the non-dogmatic SCC set initialising the model along various initial times t and the corresponding\n",
    "#   state (K, S) along the DICE trajectory\n",
    "# To run this:\n",
    "#   (5) First reinitialise for the non-dogmatic case desired, load the solution. Important to have the correct Φ loaded!\n",
    "#   (6) Then load benchmark_paths (to get the DICE trajectory)\n",
    "#   (7) Then run this block, setting calc_type below correctly\n",
    "#   (8) If computed for multiple degrees of non-dogmatism, repeat for each\n",
    "\n",
    "# This is the degree of non-dogmatism we use\n",
    "calc_type = 3;\n",
    "\n",
    "# Now compute the trajectory starting from each initial state and store the resulting SCC value, for each type\n",
    "for tt = 1:size(DICE_SCC_pairs, 3)\n",
    "    for nn = 1:size(DICE_SCC_pairs, 1)\n",
    "        zzz = get_trajectories([paths_benchmark[1][1].K_t[tt], paths_benchmark[1][1].S_t[tt]], Φ, 2, nn*ones(Integer, T), DICE_parameters, approx_parameters, functions);\n",
    "        DICE_SCC_pairs[nn, calc_type, tt] = zzz.SCC_t[1];\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Writes the resulting SCCs into a CSV file\n",
    "\n",
    "CSV.write(string(\".\\\\\", model_instance,  \"\\\\DICE_SCC.csv\"),  DataFrame(paths_benchmark[1][1].SCC_t[1:5]'), writeheader=false)\n",
    "for j = 1:5\n",
    "    CSV.write(string(\".\\\\\", model_instance,  \"\\\\DICE_SCC.csv\"),  DataFrame(DICE_SCC_pairs[:, :, j]), writeheader=false, append=true)\n",
    "end\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plots the DICE benchmark path and the non-dogmatic SCCs along the math.\n",
    "\n",
    "DICE_plot = plot(1:5, 12/44*paths_benchmark[1][1].SCC_t[1:5])\n",
    "for tt = 1:5\n",
    "    DICE_plot = scatter!(tt*ones(10).-0.1, 12/44*DICE_SCC_pairs[:, 1, tt], alpha = 1, markeralpha = 1, markerstrokealpha = 0, markerstrokecolor=base_linecolours[1:10],\n",
    "        color=base_linecolours[1:10], label=\"\", markershape = :circle, markersize = 6, xtickfont=font(7, \"Arial\"), ytickfont=font(7, \"Arial\"),\n",
    "        guidefontsize=8);\n",
    "    DICE_plot = scatter!(tt*ones(10).+0.00, 12/44*DICE_SCC_pairs[:, 2, tt], alpha = 1, markeralpha = 1, markerstrokealpha = 0, markerstrokecolor=base_linecolours[1:10],\n",
    "        color=:white, label=\"\", markershape = :circle, markersize = 6, xtickfont=font(7, \"Arial\"), ytickfont=font(7, \"Arial\"),\n",
    "        guidefontsize=8);\n",
    "    DICE_plot = scatter!(tt*ones(10).+0.10, 12/44*DICE_SCC_pairs[:, 3, tt], alpha = 1, markeralpha = 1, markerstrokealpha = 1, markerstrokecolor=base_linecolours[1:10],\n",
    "        color=base_linecolours[1:10], label=\"\", markershape = :x, markersize = 6, xtickfont=font(7, \"Arial\"), ytickfont=font(7, \"Arial\"),\n",
    "        guidefontsize=8);\n",
    "end\n",
    "\n",
    "display(DICE_plot);\n",
    "\n",
    "Plots.savefig(string(\".\\\\\", model_instance, \"\\\\DICE.pdf\"));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save the DICE-path SCC data\n",
    "\n",
    "@save string(\".\\\\\", model_instance, \"\\\\DICE_data_\", model_instance, \".jld2\") DICE_SCC_pairs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# End of task: compute DICE-path data (i.e. Figure 4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Task: plot the model outcome (in which the type never changes) for each initial type.\n",
    "#   These can then be plotted.\n",
    "#   Also save the initial SCC values and the initial consumption values.\n",
    "\n",
    "# This block sets up an array of arrays, one for each degree of dogmatism.\n",
    "# It then computes a 'trajectory' and initialises each of these arrays for the type 'trajectory'\n",
    "\n",
    "# The next line initialised the paths array the first time to start generations paths.\n",
    "# It should be commented out if you have already done some runs and want to load pre-computed data (a few blocks below)\n",
    "paths = Array{Array, 1}(undef, dogmatism_runs)\n",
    "\n",
    "zzz = get_trajectories([eff(223, 1), .5], Φ, 5, ones(Integer, T), DICE_parameters, approx_parameters, functions)\n",
    "\n",
    "paths[runnumber] = Array{typeof(zzz)}(undef, 0)\n",
    "\n",
    "for j = 1:length(βs)\n",
    "\n",
    "    zzz = get_trajectories([eff(223, 1), .5], Φ, T, j*ones(Integer, T), DICE_parameters, approx_parameters, functions)\n",
    "    push!(paths[runnumber], zzz)\n",
    "\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save the paths object for the given model\n",
    "# Note that you can then reinitialise for a different degree of non-dogmatism (runnumber), load the previously saved\n",
    "# paths variable, and run the above block again to populate a different elements of paths\n",
    "\n",
    "@save string(\".\\\\\", model_instance, \"\\\\paths_\", model_instance, \".jld2\") paths"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the paths object for the given model\n",
    "\n",
    "@load string(\".\\\\\", model_instance, \"\\\\paths_\", model_instance, \".jld2\") paths"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up a global variable \n",
    "# This again (like paths) contains the data require to plot the modal trajectory (i.e. no type changes)\n",
    "#   for each degree of non-dogmatism, and for all initial types\n",
    "\n",
    "plotdata = NamedTuple[(\n",
    "                c = zeros(T),\n",
    "                K = zeros(T),\n",
    "                S = zeros(T),\n",
    "                p = zeros(T),\n",
    "                SCC = zeros(T),\n",
    "                μ = zeros(T))\n",
    "        \n",
    "                for i in 1:dogmatism_runs]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Populate plotdata for a given run (if this needs doing)\n",
    "\n",
    "# Note that the corresponding element of the paths array must have been computed above\n",
    "\n",
    "plotdata[runnumber] = get_plotdata(paths[runnumber], T, length(βs))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Populate plotdata for all runs if desired\n",
    "\n",
    "for i = 1:dogmatism_runs\n",
    "    plotdata[i] = get_plotdata(paths[i], T, length(βs))\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save plotdata\n",
    "\n",
    "@save string(\".\\\\\", model_instance, \"\\\\plotdata_\", model_instance, \".jld2\") plotdata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Functions to plot SCC and marginal abatement cost trajectories\n",
    "\n",
    "using Plots\n",
    "using LaTeXStrings\n",
    "\n",
    "function plot_SCC(plotdata_d, plotdata_nd, T, T_zoom, linestyles, linecolours, labels)\n",
    "\n",
    "    n = size(plotdata_d.c, 1)\n",
    "        \n",
    "    plt1 = plot(1:T, (transpose(plotdata_d.SCC[:, 1:T])), title=string(model_instance, \"SCC: dogmatic\"), legend=:none, linestyle=reshape(linestyles[:,1], 1, size(linestyles, 1)), color=reshape(linecolours[:,1], 1, size(linestyles, 1)))\n",
    "    plt2 = plot(1:T, (transpose(plotdata_nd.SCC[:, 1:T])), title=\"SCC: nondogmatic\", legend=:none, linestyle=reshape(linestyles[:,2], 1, size(linestyles, 1)), color=reshape(linecolours[:,2], 1, size(linestyles, 1)))\n",
    "    plt3 = plot(1:T_zoom, (transpose(plotdata_d.SCC[:, 1:T_zoom])), ylimits=(0, 700), legend=:right, linestyle=reshape(linestyles[:,1], 1, size(linestyles, 1)), color=reshape(linecolours[:,1], 1, size(linestyles, 1)), label=reshape(labels[:,1], 1, size(linestyles, 1)))\n",
    "    plt4 = plot(1:T_zoom, (transpose(plotdata_nd.SCC[:, 1:T_zoom])), ylimits=(0, 700), legend=:right, linestyle=reshape(linestyles[:,2], 1, size(linestyles, 1)), color=reshape(linecolours[:,2], 1, size(linestyles, 1)), label=reshape(labels[:,2], 1, size(linestyles, 1)))\n",
    "\n",
    "    plot_tax = plot(plt2, plt4, plt1, plt3, layout=(2,2), size=(1000,800))\n",
    "    Plots.savefig(string(\".\\\\\", model_instance, \"\\\\SCC.pdf\"));\n",
    "    #display(plot_tax)\n",
    "    \n",
    "    return plot_tax\n",
    "end\n",
    "\n",
    "\n",
    "function plot_MAC(plotdata_d, plotdata_nd, T, T_zoom, linestyles, linecolours, labels)\n",
    "\n",
    "    n = size(plotdata_d.c, 1)\n",
    "        \n",
    "    plt1 = plot(1:T, (transpose(plotdata_d.p[:, 1:T])), title=string(model_instance, \" MAC: dogmatic\"), legend=:none, linestyle=reshape(linestyles[:,1], 1, size(linestyles, 1)), color=reshape(linecolours[:,1], 1, size(linestyles, 1)))\n",
    "    plt2 = plot(1:T, (transpose(plotdata_nd.p[:, 1:T])), title=\"MAC: nondogmatic\", legend=:none, linestyle=reshape(linestyles[:,2], 1, size(linestyles, 1)), color=reshape(linecolours[:,2], 1, size(linestyles, 1)))\n",
    "    plt3 = plot(1:T_zoom, (transpose(plotdata_d.p[:, 1:T_zoom])), ylimits=(0, 700), legend=:right, linestyle=reshape(linestyles[:,1], 1, size(linestyles, 1)), color=reshape(linecolours[:,1], 1, size(linestyles, 1)), label=reshape(labels[:,1], 1, size(linestyles, 1)))\n",
    "    plt4 = plot(1:T_zoom, (transpose(plotdata_nd.p[:, 1:T_zoom])), ylimits=(0, 700), legend=:right, linestyle=reshape(linestyles[:,2], 1, size(linestyles, 1)), color=reshape(linecolours[:,2], 1, size(linestyles, 1)), label=reshape(labels[:,2], 1, size(linestyles, 1)))\n",
    "\n",
    "    plot_tax = plot(plt2, plt4, plt1, plt3, layout=(2,2), size=(1000,800))\n",
    "    Plots.savefig(string(\".\\\\\", model_instance, \"\\\\MAC.pdf\"));\n",
    "    #display(plot_tax)\n",
    "    \n",
    "    return plot_tax\n",
    "end\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to plot capital trajectories\n",
    "\n",
    "function plot_capital(plotdata_d, plotdata_nd, T, T_zoom, linestyles, linecolours, labels)\n",
    "    \n",
    "    plt1 = plot(1:T, (transpose(plotdata_nd.K[:, 1:T])), title=string(model_instance, \" Capital: nondogmatic\"), legend=:best, linestyle=reshape(linestyles[:,2], 1, size(linestyles, 1)), color=reshape(linecolours[:,2], 1, size(linestyles, 1)), label=reshape(labels[:,2], 1, size(linestyles, 1)))\n",
    "    plt2 = plot(1:T, (transpose(plotdata_nd.K[:, 1:T])), ylimits=(0,10), title=\"Capital: nondogmatic (detail)\", legend=:none, linestyle=reshape(linestyles[:,2], 1, size(linestyles, 1)), color=reshape(linecolours[:,2], 1, size(linestyles, 1)))\n",
    "    plt3 = plot(1:T, (transpose(plotdata_d.K[:, 1:T])), title=\"Capital: dogmatic\", legend=:best, linestyle=reshape(linestyles[:,1], 1, size(linestyles, 1)), color=reshape(linecolours[:,1], 1, size(linestyles, 1)), label=reshape(labels[:,1], 1, size(linestyles, 1)))\n",
    "    #plt4 = plot(1:T, (transpose(plotdata_d.K[:, 1:T])), title=\"Capital: dogmatic (detail)\", ylimits=(0, 10), legend=:none, linestyle=[:solid :solid :solid :dash :dash :dashdot :solid :solid :dash :dash], color=[:red :green :blue :red :green :blue :black :cyan :black :cyan])\n",
    "\n",
    "    plot_K = plot(plt1, plt2, plt3, layout=(3,1), size=(1000,800))\n",
    "    Plots.savefig(string(\".\\\\\", model_instance, \"\\\\capital.pdf\"));\n",
    "    #display(plot_K)\n",
    "    \n",
    "    return plot_K\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Functions to plot abatement rate trajectories\n",
    "\n",
    "function plot_mu(plotdata_d, plotdata_nd, T, T_zoom, linestyles, linecolours, labels)\n",
    "    \n",
    "    plt1 = plot(1:T, (transpose(plotdata_nd.μ[:, 1:T])), title=string(model_instance, \" Abatement: nondogmatic\"), legend=:best, linestyle=reshape(linestyles[:,2], 1, size(linestyles, 1)), color=reshape(linecolours[:,2], 1, size(linestyles, 1)), label=reshape(labels[:,2], 1, size(linestyles, 1)))\n",
    "    #plt2 = plot(1:T, (transpose(plotdata_nd.K[:, 1:T])), ylimits=(0,10), title=\"Capital: nondogmatic (detail)\", legend=:none, linestyle=reshape(linestyles[:,2], 1, size(linestyles, 1)), color=reshape(linecolours[:,2], 1, size(linestyles, 1)))\n",
    "    plt3 = plot(1:T, (transpose(plotdata_d.μ[:, 1:T])), title=\"Abatement: dogmatic\", legend=:best, linestyle=reshape(linestyles[:,1], 1, size(linestyles, 1)), color=reshape(linecolours[:,1], 1, size(linestyles, 1)), label=reshape(labels[:,1], 1, size(linestyles, 1)))\n",
    "    #plt4 = plot(1:T, (transpose(plotdata_d.K[:, 1:T])), title=\"Capital: dogmatic (detail)\", ylimits=(0, 10), legend=:none, linestyle=[:solid :solid :solid :dash :dash :dashdot :solid :solid :dash :dash], color=[:red :green :blue :red :green :blue :black :cyan :black :cyan])\n",
    "\n",
    "    plot_K = plot(plt1, plt3, layout=(2,1), size=(1000,800))\n",
    "    Plots.savefig(string(\".\\\\\", model_instance, \"\\\\mu.pdf\"));\n",
    "    #display(plot_K)\n",
    "    \n",
    "    return plot_K\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Functions to plot consumption trajectories\n",
    "\n",
    "function plot_consumption(plotdata_d, plotdata_nd, T, T_zoom, linestyles, linecolours, labels)\n",
    "\n",
    "    plt1 = plot(1:T, (transpose(plotdata_nd.c[:, 1:T])), ylimits=(0,3e3), title=string(model_instance, \" Consumption: nondogmatic\"), legend=:none, linestyle=reshape(linestyles[:,2], 1, size(linestyles, 1)), color=reshape(linecolours[:,2], 1, size(linestyles, 1)))\n",
    "    plt2 = plot(1:T_zoom, (transpose(plotdata_nd.c[:, 1:T_zoom])), ylimits=(20,120), title=\"(detail)\", legend=:best, linestyle=reshape(linestyles[:,2], 1, size(linestyles, 1)), color=reshape(linecolours[:,2], 1, size(linestyles, 1)), label=reshape(labels[:,2], 1, size(linestyles, 1)))\n",
    "    #plt12 = plot(1:T_zoom, (transpose(plotdata_nd.c[:, 1:T_zoom])), legend = :right, label=reshape(labels[:,2], 1, 10))\n",
    "    plt3 = plot(1:T, (transpose(plotdata_d.c[:, 1:T])), ylimits=(0,3e3), title=\"Consumption: dogmatic\", legend=:none, linestyle=reshape(linestyles[:,1], 1, size(linestyles, 1)), color=reshape(linecolours[:,1], 1, size(linestyles, 1)))\n",
    "    plt4 = plot(1:T_zoom, (transpose(plotdata_d.c[:, 1:T_zoom])), ylimits=(20, 120), title=\"(detail)\", legend=:best, linestyle=reshape(linestyles[:,1], 1, size(linestyles, 1)), color=reshape(linecolours[:,1], 1, size(linestyles, 1)), label=reshape(labels[:,1], 1, size(linestyles, 1)))\n",
    "    #plt34 = plot(1:T_zoom, (transpose(plotdata_d.c[:, 1:T_zoom])), legend = :right, label=reshape(labels[:,1], 1, 10))\n",
    "\n",
    "    plot_K = plot(plt1, plt2, plt3, plt4, layout=(2,2), size=(1000,600))\n",
    "    Plots.savefig(string(\".\\\\\", model_instance, \"\\\\consumption.pdf\"));\n",
    "    #display(plot_K)\n",
    "    \n",
    "    return plot_K\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load pre-computed plotdata if you want to\n",
    "\n",
    "@load string(\".\\\\\", model_instance, \"\\\\plotdata_\", model_instance, \".jld2\") plotdata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up some graphical parameters\n",
    "\n",
    "using Colors;\n",
    "\n",
    "base_linestyles = [:solid :solid :solid :solid :solid :dash :dash :dash :dash :dash]\n",
    "base_linecolours = [RGB(0.0,0.0,1.0) RGB(1.0,0.0,0.0) RGB(0.0,1.0,0.0) RGB(0.0,0.0,0.1724) RGB(1.0,0.1034,0.7241) RGB(1.0,0.8276,0.0) RGB(0.0,0.3448,0.0) RGB(0.5172,0.5172,1.0) RGB(0.6207,0.3103,0.2759) RGB(0.0,1.0,0.7586) :red :green :blue :black :cyan]\n",
    "\n",
    "years = (2015:10:2565)\n",
    "\n",
    "linestyles = Array{Symbol, 2}(undef, length(modelrun), length(βs))\n",
    "linecolours = Array{RGB{Float64}, 2}(undef, length(modelrun), length(βs))\n",
    "labels = Array{String, 2}(undef, length(modelrun), length(βs))\n",
    "\n",
    "for i in 1:length(modelrun)\n",
    "    for (index, position) in enumerate(sortperm(modelrun[i].ηs))\n",
    "        #global type_indices = sortperm(modelrun[i].ηs)\n",
    "        linestyles[i, position] = base_linestyles[index]\n",
    "        linecolours[i, position] = base_linecolours[index]\n",
    "        labels[i, position] = \"β = $(modelrun[i].βs[position]), η = $(modelrun[i].ηs[position])\"\n",
    "    end\n",
    "end\n",
    "\n",
    "labels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot capital\n",
    "\n",
    "pyplot()\n",
    "plot_capital(plotdata[1], plotdata[4], 50, 6, [linestyles[1, :] linestyles[4,:]], [linecolours[1, :] linecolours[4,:]], [labels[1,:] labels[4,:]]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot consumption\n",
    "\n",
    "pyplot()\n",
    "plot_consumption(plotdata[1], plotdata[4], T, 6, [linestyles[1, :] linestyles[4,:]], [linecolours[1, :] linecolours[4,:]], [labels[1,:] labels[4,:]]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot SCC\n",
    "\n",
    "pyplot()\n",
    "plot_SCC(plotdata[1], plotdata[4], T, 10, [linestyles[1, :] linestyles[4,:]], [linecolours[1, :] linecolours[4,:]], [labels[1,:] labels[4,:]]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot marginal abatement cost\n",
    "\n",
    "pyplot()\n",
    "plot_MAC(plotdata[1], plotdata[2], T, 10, [linestyles[1, :] linestyles[4,:]], [linecolours[1, :] linecolours[4,:]], [labels[1,:] labels[4,:]]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot abatement rates\n",
    "\n",
    "pyplot()\n",
    "plot_mu(plotdata[1], plotdata[1], T, 6, [linestyles[1, :] linestyles[4,:]], [linecolours[1, :] linecolours[4,:]], [labels[1,:] labels[4,:]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute the initial (period-0) SCC for all the types and save\n",
    "\n",
    "using CSV, DataFrames;\n",
    "\n",
    "initial_SCC = zeros(dogmatism_runs, length(βs))\n",
    "\n",
    "for i in 1:dogmatism_runs\n",
    "    if size(plotdata[i].SCC, 2) > 1\n",
    "        initial_SCC[i, :] = plotdata[i].SCC[:,1];\n",
    "    end\n",
    "end\n",
    "\n",
    "CSV.write(string(\".\\\\\", model_instance,  \"\\\\initial_SCC.csv\"),  DataFrame(initial_SCC), writeheader=false)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute the initial (period-0) consumption rates for all the types and save\n",
    "\n",
    "using CSV, DataFrames;\n",
    "\n",
    "initial_c = zeros(dogmatism_runs, length(βs))\n",
    "\n",
    "for i in 1:dogmatism_runs\n",
    "    if size(plotdata[i].c, 2) > 1\n",
    "        initial_c[i, :] = plotdata[i].c[:,1];\n",
    "    end\n",
    "end\n",
    "\n",
    "CSV.write(string(\".\\\\\", model_instance,  \"\\\\initial_c.csv\"),  DataFrame(initial_c), writeheader=false)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# End of task: plot modal runs, save initial SCC and consumption values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Task: compute Monte Carlo runs\n",
    "# Here we initialise many realisations of the type transition process and compute the actual trajectory of\n",
    "#   the economy for each realisation.\n",
    "\n",
    "# Calculate the cumulative weight transition matrix\n",
    "\n",
    "cumulative_weights = [ zeros(agent_type_n, agent_type_n) for i in 1:dogmatism_runs ];\n",
    "\n",
    "for i in 1:agent_type_n\n",
    "    cumulative_weights[runnumber][:, i] = sum(weights[:, 1:i], dims=2);\n",
    "end\n",
    "\n",
    "cumulative_weights[runnumber]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up fundamental fan chart / Monte Carlo parameters\n",
    "\n",
    "fan_runs = 10;\n",
    "fan_initvalue = [1 2 3 4 5 6 7 8 9 10];\n",
    "n_draws = 1;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load pre-generated data for fan chart / Monter Carlo run\n",
    "# Note: generating the fan_paths variable is quite computationally intensive. Recommended to load it if possible. But the below code will\n",
    "#   also generate it from scratch if desired.\n",
    "\n",
    "@load string(\".\\\\\", model_instance, \"\\\\fandata_\", model_instance, \".jld2\") fan_runs fan_initvalue n_draws draws typechain fan_paths"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate the Markov chain of types *if* starting from scratch...\n",
    "#   If the fan/MC data have already been generated, they can be loaded using a block below.\n",
    "\n",
    "draws = [ zeros(n_draws, T-1) for i in 1:fan_runs ]\n",
    "typechain = [ zeros(Int8, n_draws, T) for i in 1:fan_runs]\n",
    "for i in 1:fan_runs\n",
    "    draws[i] = rand(Float64, n_draws, T-1)\n",
    "    typechain[i][:, 1] = fan_initvalue[i]*ones(n_draws);\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate realisations of the type path\n",
    "\n",
    "for j in 1:fan_runs\n",
    "    for t in 1:T-1\n",
    "        for i in 1:n_draws\n",
    "            typechain[j][i, t+1] = min(10, searchsortedfirst(cumulative_weights[runnumber][typechain[j][i, t],:], draws[j][i, t]));\n",
    "        end\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialise the fan paths\n",
    "\n",
    "fan_paths = Array{NamedTuple, 2}(undef, fan_runs, n_draws)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now calculate all the paths\n",
    "# Note: At least my implementation had a memory leak or something, which meant computation becomes exceedingly slow after a few hundred\n",
    "#   runs (at least if data for several types has already been computed and is being stored).  In case of slow execution, may be advisable\n",
    "#   to interrupt, save the computed data, restart and reinitialise everything, and then start computing again from the interruption on.\n",
    "\n",
    "for fan_runnumber in 1:10\n",
    "    \n",
    "    for i in 1:n_draws\n",
    "        println(\"Draw $i\")\n",
    "        zzz = get_trajectories([eff(223, 1), .5], Φ, 50, typechain[fan_runnumber][i, :], DICE_parameters, approx_parameters, functions)\n",
    "        fan_paths[fan_runnumber, i] = zzz;\n",
    "    end\n",
    "end\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save the fan / Monte Carlo data\n",
    "\n",
    "@save string(\".\\\\\", model_instance, \"\\\\fandata_\", model_instance, \".jld2\") fan_runs fan_initvalue n_draws draws typechain fan_paths"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialise the 'spaghetti' plot array which contains the data for each Monte Carlo run\n",
    "# This must be generated even if the fan_paths above has been loaded.\n",
    "# The spaghetti data is crucial to compute the distributions underlying Figure 3.\n",
    "\n",
    "spaghetti = Array{NamedTuple, 2}(undef, fan_runs, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now get the actual plot data\n",
    "\n",
    "for i = 1:10\n",
    "    spaghetti[i] = get_plotdata(fan_paths[i, :], T, n_draws)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot histograms of capital distributions in a given time period.\n",
    "#   This is ad hoc coding for internal analysis purposes but retained so it can be seen how\n",
    "#   the process works\n",
    "\n",
    "distr_plot_1 = histogram(spaghetti_3.K[:, 30], bins=0:5:50, label=\"K_3 (30)\");\n",
    "distr_plot_2 = histogram(spaghetti_4.K[:, 30], bins=0:5:50, label=\"K_4 (30)\");\n",
    "distr_plot_3 = histogram(spaghetti_3.S[:, 30], bins=0.6:.1:1.6, label=\"S_3 (30)\");\n",
    "distr_plot_4 = histogram(spaghetti_4.S[:, 30], bins=0.6:.1:1.6, label=\"S_4 (30)\");\n",
    "plot(distr_plot_1, distr_plot_2, distr_plot_3, distr_plot_4)\n",
    "Plots.savefig(string(\".\\\\\", model_instance, \"\\\\ergodicity_34.pdf\"));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute some means and standard deviations from the Monte Carlo exercise.\n",
    "\n",
    "spaghetti_stats = Array{NamedTuple, 1}(undef, fan_runs);\n",
    "\n",
    "for i = 1:fan_runs\n",
    "    spaghetti_stats[i] = (g_c_mean = Array{Float64, 1}(undef, T), g_c_sd = Array{Float64, 1}(undef, T), S_mean = Array{Float64, 1}(undef, T), K_mean = Array{Float64, 1}(undef, T), r_mean = Array{Float64, 1}(undef, T), r_sd = Array{Float64, 1}(undef, T));\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute the consumption discount rates (and components of this).\n",
    "\n",
    "r_chi = [zeros(n_draws, T) for _ in 1:10]\n",
    "delta_chi = [zeros(n_draws, T) for _ in 1:10]\n",
    "r_chi = [zeros(n_draws, T) for _ in 1:10]\n",
    "rho_chi = [zeros(n_draws, T) for _ in 1:10]\n",
    "G_chi = [zeros(n_draws, T) for _ in 1:10]\n",
    "P_chi = [zeros(n_draws, T) for _ in 1:10]\n",
    "MU_chi = [zeros(n_draws, T) for _ in 1:10]\n",
    "MD_chi = [zeros(n_draws, T) for _ in 1:10]\n",
    "\n",
    "for i in 1:fan_runs\n",
    "    delta_chi[i][:,1] = ones(n_draws, 1);\n",
    "    for j in 1:n_draws\n",
    "        for k in 1:T\n",
    "            if k > 1\n",
    "                delta_chi[i][j,k] = delta_chi[i][j,k-1] .* βs[typechain[i][j,k-1]];\n",
    "            end\n",
    "            r_chi[i][j,k] = (L_seq[k]*delta_chi[i][j,k]*du(spaghetti[i].c[j,k]/L_seq[k], ηs[typechain[i][j,k]])/(L_seq[1]*du(spaghetti[i].c[j,1]/L_seq[1], ηs[typechain[i][j,1]])))^(-1/((k-1)*period_length)) - 1;\n",
    "            rho_chi[i][j,k] = (delta_chi[i][j,k])^(-1/((k-1)*period_length)) - 1;\n",
    "            G_chi[i][j,k] = (du(spaghetti[i].c[j,k]/L_seq[k], ηs[typechain[i][j,k]])/(du(spaghetti[i].c[j,1]/L_seq[1], ηs[typechain[i][j,1]])))^(-1/((k-1)*period_length)) - 1;\n",
    "            P_chi[i][j,k] = (L_seq[k]/L_seq[1])^(1/((k-1)*period_length)) - 1;\n",
    "            MU_chi[i][j,k] = L_seq[k].*du(spaghetti[i].c[j,k]/L_seq[k], ηs[typechain[i][j,k]]);\n",
    "            MD_chi[i][j,k] = dDamage(spaghetti[i].S[j,k]) .* F(A_seq[k], spaghetti[i].K[j,k], L_seq[k]) * 44/12;\n",
    "        end\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute the mean carbon tax (proxied by the marginal damage)\n",
    "\n",
    "mean_tax = zeros(agent_type_n, T)\n",
    "for i = 1:agent_type_n\n",
    "    for k = 1:T\n",
    "        mean_tax[i,k] = mean(MD_chi[i][:,k]);\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute the means of the discount rates and insert into the Monte Carlo statistics array\n",
    "\n",
    "using Statistics;\n",
    "\n",
    "t = 31;\n",
    "\n",
    "for i = 1:10\n",
    "    MC_mean = mean((spaghetti[i].c[:, t]./spaghetti[i].c[:, t-1]).^(1/period_length) - ones(1000,1));\n",
    "    MC_sd = std((spaghetti[i].c[:, t]./spaghetti[i].c[:, t-1]).^(1/period_length) - ones(1000,1));\n",
    "    MC_S_mean = mean(spaghetti[i].S[:, t])\n",
    "    MC_K_mean = mean(spaghetti[i].K[:, t])\n",
    "    MC_discount_mean = mean(r_chi[i][:, t]);\n",
    "    MC_discount_sd = std(r_chi[i][:, t]);\n",
    "    spaghetti_stats[i].g_c_mean[t] = MC_mean;\n",
    "    spaghetti_stats[i].g_c_sd[t] = MC_sd;\n",
    "    spaghetti_stats[i].S_mean[t] = MC_S_mean;    \n",
    "    spaghetti_stats[i].K_mean[t] = MC_K_mean;    \n",
    "    spaghetti_stats[i].r_mean[t] = MC_discount_mean;\n",
    "    spaghetti_stats[i].r_sd[t] = MC_discount_sd;\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# For the same t as used above, print the statistics by type\n",
    "\n",
    "using Printf;\n",
    "\n",
    "println(\"Absolute growth rates:\")\n",
    "print(\"  Mean:\");\n",
    "for i = 1:10\n",
    "    print(@sprintf(\" %.4f\", spaghetti_stats[i].g_c_mean[t]));\n",
    "end\n",
    "println();\n",
    "print(\"  Stdev:\")\n",
    "for i = 1:10\n",
    "    print(@sprintf(\" %.4f\", spaghetti_stats[i].g_c_sd[t]));\n",
    "end\n",
    "println();\n",
    "\n",
    "print(\"CO2 stock:\")\n",
    "for i = 1:10\n",
    "    print(@sprintf(\" %.4f\", spaghetti_stats[i].S_mean[t]));\n",
    "end\n",
    "println();\n",
    "\n",
    "print(\"Capital stock:\")\n",
    "for i = 1:10\n",
    "    print(@sprintf(\" %.4f\", spaghetti_stats[i].K_mean[t]));\n",
    "end\n",
    "println();\n",
    "\n",
    "println(\"Consumption discount rates:\")\n",
    "print(\"  Mean:\");\n",
    "for i = 1:10\n",
    "    print(@sprintf(\" %.4f\", spaghetti_stats[i].r_mean[t]));\n",
    "end\n",
    "println();\n",
    "print(\"  Stdev:\")\n",
    "for i = 1:10\n",
    "    print(@sprintf(\" %.4f\", spaghetti_stats[i].r_sd[t]));\n",
    "end\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot histograms of consumption growth rates in a given period\n",
    "\n",
    "distr_plot = Array{Plots.Plot,1}(undef,20)\n",
    "plot_types = [1 2 3 4 5 6 7 8  9 10];\n",
    "t = 40;\n",
    "\n",
    "for i = 1:length(plot_types)\n",
    "    distr_plot[i] = histogram((spaghetti[plot_types[i]].c[:, t]./spaghetti[plot_types[i]].c[:, t-1]).^(1/period_length) - ones(1000,1), bins=-.04:.002:.04, label=\"g_c_$(plot_types[i]) ($t)\", ylims=(0,400));\n",
    "end\n",
    "\n",
    "plot(distr_plot[1:length(plot_types)]...)\n",
    "Plots.savefig(string(\".\\\\\", model_instance, \"\\\\ergodicity_cons.pdf\"));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot histograms of consumption discount rates in a given period\n",
    "\n",
    "distr_plot = Array{Plots.Plot,1}(undef,20)\n",
    "plot_types = [1 2 3 4 5 6 7 8  9 10];\n",
    "t = 40;\n",
    "\n",
    "for i = 1:length(plot_types)\n",
    "    distr_plot[i] = histogram(r_chi[i][:, t], ylims=(0,500), bins=0:.025:1, legend=false);\n",
    "end\n",
    "\n",
    "plot(distr_plot[1:length(plot_types)]...)\n",
    "Plots.savefig(string(\".\\\\\", model_instance, \"\\\\ergodicity_discount.pdf\"));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot histograms of marginal damages in a given period\n",
    "\n",
    "distr_plot = Array{Plots.Plot,1}(undef,20)\n",
    "plot_types = [1 2 3 4 5 6 7 8  9 10];\n",
    "t = 31;\n",
    "\n",
    "for i = 1:length(plot_types)\n",
    "    distr_plot[i] = Plots.histogram(MD_chi[i][:, t], ylims=(0,500), legend=false);\n",
    "end\n",
    "\n",
    "Plots.plot(distr_plot[1:length(plot_types)]...)\n",
    "Plots.savefig(string(\".\\\\\", model_instance, \"\\\\ergodicity_payoff_histogram.pdf\"));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the discount rate term structure figure, as well as the components\n",
    "\n",
    "using Colors;\n",
    "using GRUtils;\n",
    "using Measures;\n",
    "using LaTeXStrings;\n",
    "using PyPlot;\n",
    "using Statistics;\n",
    "\n",
    "pyplot();\n",
    "Plots.PyPlotBackend();\n",
    "\n",
    "base_linecolours = [RGB(0.0,0.0,1.0) RGB(1.0,0.0,0.0) RGB(0.0,1.0,0.0) RGB(0.0,0.0,0.1724) RGB(1.0,0.1034,0.7241) RGB(1.0,0.8276,0.0) RGB(0.0,0.3448,0.0) RGB(0.5172,0.5172,1.0) RGB(0.6207,0.3103,0.2759) RGB(0.0,1.0,0.7586) :red :green :blue :black :cyan]\n",
    "\n",
    "plot_times = [2 6 11 16 21 26 31];\n",
    "#plot_times = [2 3 4 5 6 7 8];\n",
    "\n",
    "r_plot = Plots.plot(xlims=(0, 5*(length(plot_times)+1)-6.5), ylims=(0, 10), legend=:topright, legendfontsize=6, fontfamily=\"Helvetica\", title=\"Consumption discount rate (\\$r(t|χ)\\$)\", ylabel=\"%/yr\", xlabel=\"Maturity (yrs)\", grid=false, xtickfont=font(9, \"Helvetica\"), ytickfont=font(9, \"Helvetica\"),\n",
    "        guidefontsize=9, titlefontfamily=\"Helvetica\", titlefontsize=9, guidefontfamily=\"Helvetica\", foreground_color_legend = :white, background_color_legend = :transparent);\n",
    "# For the composite plot: bottom_margin=10mm\n",
    "\n",
    "xdata = zeros(length(plot_times));\n",
    "ydata = zeros(length(plot_times));\n",
    "errordata = zeros(length(plot_times), 5);\n",
    "\n",
    "rho_ydata = zeros(length(plot_times));\n",
    "rho_errordata = zeros(length(plot_times), 5);\n",
    "G_ydata = zeros(length(plot_times));\n",
    "G_errordata = zeros(length(plot_times), 5);\n",
    "P_ydata = zeros(length(plot_times));\n",
    "P_errordata = zeros(length(plot_times), 5);\n",
    "MU_ydata = zeros(length(plot_times));\n",
    "MU_errordata = zeros(length(plot_times), 5);\n",
    "MD_ydata = zeros(length(plot_times));\n",
    "MD_errordata = zeros(length(plot_times), 5);\n",
    "\n",
    "\n",
    "rho_plot = Plots.plot(xlims=(0,5*(length(plot_times)+1)-6.5), ylims=(0, 8), legend=:none, legendfontsize=6, foreground_color_legend = :white, background_color_legend = :transparent, title=\"Pure rate of time preference (\\$\\\\bar{\\\\rho}_{t|\\\\chi}\\$)\", ylabel=\"%/yr\", xlabel=\"Maturity (yrs)\", grid=false, xtickfont=font(9, \"Helvetica\"), ytickfont=font(9, \"Helvetica\"),\n",
    "        guidefontsize=9, titlefontfamily=\"Helvetica\", titlefontsize=9, guidefontfamily=\"Helvetica\");\n",
    "#rho_plot = plot(xlims=(0,5*(length(plot_times)+1)-6.5), ylims=(0, 8), legend=:none, title=\"Pure rate of time preference\", grid=false, xtickfont=font(9, \"Helvetica\"), ytickfont=font(9, \"Helvetica\"),\n",
    "#        guidefontsize=9, titlefontfamily=\"Helvetica\", titlefontsize=9, guidefontfamily=\"Helvetica\");\n",
    "G_plot = Plots.plot(xlims=(0,5*(length(plot_times)+1)-6.5), ylims=(0, 8), legend=:none, legendfontsize=6, title=\"Consumption smoothing (\\$G_{t|\\\\chi}\\$)\", ylabel=\"%/yr\", xlabel=\"Maturity (yrs)\", foreground_color_legend = :white, background_color_legend = :transparent, grid=false, xtickfont=font(9, \"Helvetica\"), ytickfont=font(9, \"Helvetica\"),\n",
    "        guidefontsize=9, titlefontfamily=\"Helvetica\", titlefontsize=9, guidefontfamily=\"Helvetica\");\n",
    "P_plot = Plots.plot(xlims=(0,5*(length(plot_times)+1)-6.5), ylims=(0, 10), legend=:none, ylabel=\"population growth\", xlabel=\"Maturity\", grid=false, xtickfont=font(9, \"Helvetica\"), ytickfont=font(9, \"Helvetica\"),\n",
    "        guidefontsize=9, titlefontfamily=\"Helvetica\", titlefontsize=9, guidefontfamily=\"Helvetica\");\n",
    "MU_plot = Plots.plot(xlims=(0,5*(length(plot_times)+1)-6.5), legend=:none, ylabel=\"Payoff flow\", xlabel=\"Maturity\", grid=false, xtickfont=font(9, \"Helvetica\"), ytickfont=font(9, \"Helvetica\"),\n",
    "        guidefontsize=9, titlefontfamily=\"Helvetica\", titlefontsize=9, guidefontfamily=\"Helvetica\");\n",
    "MD_plot = Plots.plot(xlims=(0,5*(length(plot_times)+1)-6.5), ylims=(0, 45), legend=:topleft, legendfontsize=7, fontfamily=\"Helvetica\", ylabel=\"Marginal damage (2010\\$/tCO₂)\", xlabel=\"Maturity (yrs)\", grid=false, xtickfont=font(9, \"Helvetica\"), ytickfont=font(9, \"Helvetica\"),\n",
    "        guidefontsize=9, titlefontfamily=\"Helvetica\", titlefontsize=9, guidefontfamily=\"Helvetica\", foreground_color_legend = :white, background_color_legend = :transparent);\n",
    "\n",
    "for i = 1:fan_runs\n",
    "    for t in 1:length(plot_times)\n",
    "        #xdata[t] = ((plot_times[t] - 1) - 1.2 + (i-1)/(fan_runs-1) * 2.4)*10;\n",
    "        xdata[t] = (t*5 - 1.5 + (i-1)/(fan_runs-1) * 3.0);\n",
    "        if t > 1\n",
    "            xdata[t] = xdata[t] - 3.5;\n",
    "        else\n",
    "            xdata[t] = xdata[t] - 2.5;\n",
    "        end\n",
    "        ydata[t] = mean(r_chi[i][:, plot_times[t]])*100;\n",
    "        errordata[t, :] = quantile(r_chi[i][:, plot_times[t]], [0.05 .25 .5 .75 .95])*100;\n",
    "        if errordata[t, 2] == errordata[t, 4]\n",
    "            errordata[t,2] = NaN;\n",
    "            errordata[t,4] = NaN;\n",
    "        end\n",
    "        rho_ydata[t] = mean(rho_chi[i][:, plot_times[t]])*100;\n",
    "        rho_errordata[t, :] = quantile(rho_chi[i][:, plot_times[t]], [0.05 .25 .5 .75 .95])*100;\n",
    "        G_ydata[t] = mean(G_chi[i][:, plot_times[t]])*100;\n",
    "        G_errordata[t, :] = quantile(G_chi[i][:, plot_times[t]], [0.05 .25 .5 .75 .95])*100;\n",
    "        if G_errordata[t, 2] == G_errordata[t, 4]\n",
    "            G_errordata[t,2] = NaN;\n",
    "            G_errordata[t,4] = NaN;\n",
    "        end\n",
    "        P_ydata[t] = mean(P_chi[i][:, plot_times[t]])*100;\n",
    "        P_errordata[t, :] = quantile(P_chi[i][:, plot_times[t]], [.05 .25 .5 .75 .95])*100;\n",
    "        #MU_ydata[t] = mean(MU_chi[i][:, plot_times[t]]);\n",
    "        #MU_errordata[t, :] = quantile(MU_chi[i][:, plot_times[t]], [.05 .25 .5 .75 .95]);\n",
    "        MU_ydata[t] = mean(log.(MU_chi[i][:, plot_times[t]]));\n",
    "        MU_errordata[t, :] = quantile(log.(MU_chi[i][:, plot_times[t]]), [.05 .25 .5 .75 .95]);\n",
    "        MD_ydata[t] = mean(MD_chi[i][:, plot_times[t]]);\n",
    "        MD_errordata[t, :] = quantile(MD_chi[i][:, plot_times[t]], [.05 .25 .5 .75 .95]);\n",
    "    end\n",
    "    Plots.scatter!(r_plot, xdata, (errordata[:,1] + errordata[:, 5])/2, yerror = (errordata[:,5] - errordata[:,1])/2, markershape=:none, markersize=0, linecolor=base_linecolours[i],  label=\"\");\n",
    "    Plots.scatter!(r_plot, xdata, (errordata[:,2] + errordata[:, 4])/2, yerror = (errordata[:,4] - errordata[:,2])/2, markershape=:none, markersize=0, linecolor=base_linecolours[i],  label=\"\");\n",
    "    #Plots.scatter!(r_plot, xdata, (errordata[:,1] + errordata[:, 5])/2, yerror = (errordata[:,5] - errordata[:,1])/2, markershape=:none, markersize=0, linecolor=base_linecolours[i],  markerstrokewidth=1.25, label=\"\");\n",
    "    #Plots.scatter!(r_plot, xdata, (errordata[:,2] + errordata[:, 4])/2, yerror = (errordata[:,4] - errordata[:,2])/2, markershape=:none, markersize=0, linecolor=base_linecolours[i],  markerstrokewidth=2.5, label=\"\");\n",
    "    ##scatter!(r_plot, xdata, errordata[:,2], ylabel=\" \", markershape=:hline, markersize=3, linecolor=base_linecolours[i], markercolor=base_linecolours[i], markerstrokecolor=base_linecolours[i], markerstrokewidth=1.0, label=\"\");\n",
    "    ##scatter!(r_plot, xdata, errordata[:,4], ylabel=\" \", markershape=:hline, markersize=3, linecolor=base_linecolours[i],  markercolor=base_linecolours[i], markerstrokecolor=base_linecolours[i], markerstrokewidth=1.0, label=\"\");\n",
    "    ##scatter!(r_plot, xdata, errordata[:,3], ylabel=\" \", markershape=:diamond, markersize=3, linecolor=base_linecolours[i],  markercolor=base_linecolours[i], markerstrokecolor=base_linecolours[i], markerstrokewidth=1.0, label=\"\");\n",
    "    Plots.scatter!(r_plot, xdata, ydata, markersize=3.5, markercolor=base_linecolours[i], markerstrokewidth=0, markerstrokecolor=base_linecolours[i],\n",
    "        #label=L\"\\rho\\mathrm{ = $(round(100*(modelrun[runnumber].βs[i].^(-.1) - 1), digits=1))%, }\\eta\\mathrm{ = $(round(modelrun[runnumber].ηs[i], digits=1))}\");\n",
    "        label=\"\\$\\\\rho\\$  =  $(round(100*(modelrun[runnumber].βs[i].^(-.1) - 1), digits=1))%,  \\$\\\\eta\\$  =  $(round(modelrun[runnumber].ηs[i], digits=1))\");\n",
    "        #label=\"\\$\\\\rho \\\\textrm{ = }\\$\");\n",
    "        #);\n",
    "    #legend(\"a\", \"b\", \"a\", \"b\", \"a\", \"b\", \"a\", \"b\", \"a\", \"b\", \"a\", \"b\", \"a\", \"b\", \"a\", \"b\", \"a\", \"b\", \"a\", \"b\", maxrows=2);\n",
    "\n",
    "    Plots.xticks!(vec([2.5 6.5 11.5 16.5 21.5 26.5 31.5]), string.(vec(10 .* (plot_times .- 1))));\n",
    "    Plots.scatter!(rho_plot, xdata, (rho_errordata[:,1] + rho_errordata[:, 5])/2, yerror = (rho_errordata[:,5] - rho_errordata[:,1])/2, markershape=:none, markersize=0, linecolor=base_linecolours[i],  markerstrokewidth=.75, label=\"\");\n",
    "    Plots.scatter!(rho_plot, xdata, (rho_errordata[:,2] + rho_errordata[:, 4])/2, yerror = (rho_errordata[:,4] - rho_errordata[:,2])/2, markershape=:none, markersize=0, linecolor=base_linecolours[i],  markerstrokewidth=1.25, label=\"\");\n",
    "\n",
    "    #Plots.scatter!(rho_plot, xdata, rho_ydata, markersize=2.5, markercolor=base_linecolours[i], markerstrokewidth=0, markerstrokecolor=base_linecolours[i]);\n",
    "    Plots.scatter!(rho_plot, xdata, rho_ydata, markersize=2.5, markercolor=base_linecolours[i], markerstrokewidth=0, markerstrokecolor=base_linecolours[i],\n",
    "        #label=L\"\\rho\\mathrm{ = $(round(100*(modelrun[runnumber].βs[i].^(-.1) - 1), digits=1))%, }\\eta\\mathrm{ = $(round(modelrun[runnumber].ηs[i], digits=1))}\");\n",
    "        label=\"\\$\\\\rho\\$  =  $(round(100*(modelrun[runnumber].βs[i].^(-.1) - 1), digits=1))%,  \\$\\\\eta\\$  =  $(round(modelrun[runnumber].ηs[i], digits=1))\");\n",
    "    Plots.xticks!(vec([2.5 6.5 11.5 16.5 21.5 26.5 31.5]), string.(vec(10 .* (plot_times .- 1))));\n",
    "    Plots.scatter!(G_plot, xdata, (G_errordata[:,1] + G_errordata[:, 5])/2, yerror = (G_errordata[:,5] - G_errordata[:,1])/2, markershape=:none, markersize=0, linecolor=base_linecolours[i],  markerstrokewidth=.75, label=\"\");\n",
    "    Plots.scatter!(G_plot, xdata, (G_errordata[:,2] + G_errordata[:, 4])/2, yerror = (G_errordata[:,4] - G_errordata[:,2])/2, markershape=:none, markersize=0, linecolor=base_linecolours[i],  markerstrokewidth=1.25, label=\"\");\n",
    "    #Plots.scatter!(G_plot, xdata, G_ydata, markersize=2.5, markercolor=base_linecolours[i], markerstrokewidth=0, markerstrokecolor=base_linecolours[i]);\n",
    "    Plots.scatter!(G_plot, xdata, G_ydata, markersize=2.5, markercolor=base_linecolours[i], markerstrokewidth=0, markerstrokecolor=base_linecolours[i],\n",
    "        #label=L\"\\rho\\mathrm{ = $(round(100*(modelrun[runnumber].βs[i].^(-.1) - 1), digits=1))%, }\\eta\\mathrm{ = $(round(modelrun[runnumber].ηs[i], digits=1))}\");\n",
    "        label=\"\\$\\\\rho\\$  =  $(round(100*(modelrun[runnumber].βs[i].^(-.1) - 1), digits=1))%,  \\$\\\\eta\\$  =  $(round(modelrun[runnumber].ηs[i], digits=1))\");\n",
    "    Plots.xticks!(vec([2.5 6.5 11.5 16.5 21.5 26.5 31.5]), string.(vec(10 .* (plot_times .- 1))));\n",
    "    Plots.scatter!(P_plot, xdata, (P_errordata[:,1] + P_errordata[:, 2])/2, yerror = (P_errordata[:,2] - P_errordata[:,1])/2, markershape=:none, markersize=0, linecolor=base_linecolours[i],  markerstrokewidth=1.0);\n",
    "    Plots.scatter!(P_plot, xdata, P_ydata, markersize=2.5, markercolor=base_linecolours[i], markerstrokewidth=0, markerstrokecolor=base_linecolours[i]);\n",
    "    Plots.xticks!(vec([2.5 6.5 11.5 16.5 21.5 26.5 31.5]), string.(vec(10 .* (plot_times .- 1))));\n",
    "    Plots.scatter!(MU_plot, xdata, (MU_errordata[:,1] + MU_errordata[:, 5])/2, yerror = (MU_errordata[:,5] - MU_errordata[:,1])/2, markershape=:none, markersize=0, linecolor=base_linecolours[i],  markerstrokewidth=1.0);\n",
    "    Plots.scatter!(MU_plot, xdata, (MU_errordata[:,2] + MU_errordata[:, 4])/2, yerror = (MU_errordata[:,4] - MU_errordata[:,2])/2, markershape=:none, markersize=0, linecolor=base_linecolours[i],  markerstrokewidth=1.25, label=\"\");\n",
    "    Plots.scatter!(MU_plot, xdata, MU_ydata, markersize=2.5, markercolor=base_linecolours[i], markerstrokewidth=0, markerstrokecolor=base_linecolours[i]);\n",
    "    Plots.xticks!(vec([2.5 6.5 11.5 16.5 21.5 26.5 31.5]), string.(vec(10 .* (plot_times .- 1))));\n",
    "    Plots.scatter!(MD_plot, xdata, (MD_errordata[:,1] + MD_errordata[:, 5])/2, yerror = (MD_errordata[:,5] - MD_errordata[:,1])/2, markershape=:none, markersize=0, linecolor=base_linecolours[i],  markerstrokewidth=1.25, label=\"\");\n",
    "    Plots.scatter!(MD_plot, xdata, (MD_errordata[:,2] + MD_errordata[:, 4])/2, yerror = (MD_errordata[:,4] - MD_errordata[:,2])/2, markershape=:none, markersize=0, linecolor=base_linecolours[i],  markerstrokewidth=2.5, label=\"\");\n",
    "    Plots.scatter!(MD_plot, xdata, MD_ydata, markersize=2.5, markercolor=base_linecolours[i], markerstrokewidth=0, markerstrokecolor=base_linecolours[i], label=\"\\$\\\\rho\\$  =  $(round(100*(modelrun[runnumber].βs[i].^(-.1) - 1), digits=1))%,  \\$\\\\eta\\$  =  $(round(modelrun[runnumber].ηs[i], digits=1))\");\n",
    "    Plots.xticks!(vec([2.5 6.5 11.5 16.5 21.5 26.5 31.5]), string.(vec(10 .* (plot_times .- 1))));\n",
    "end\n",
    "\n",
    "ergodicity_layout = @layout [ a{.6h}; [b c] ]\n",
    "\n",
    "erg_plot = Plots.plot(r_plot, rho_plot, G_plot, layout = ergodicity_layout, thickness_scaling=.8, size=(600,700));\n",
    "display(erg_plot)\n",
    "Plots.savefig(erg_plot, string(\".\\\\\", model_instance, \"\\\\ergodicity_convergence_all.eps\"));\n",
    "\n",
    "ergodicity_layout_alt = @layout [ a{.7w} [b; c] ]\n",
    "\n",
    "erg_plot_alt = Plots.plot(r_plot, rho_plot, G_plot, layout = ergodicity_layout_alt, thickness_scaling=.8);\n",
    "display(erg_plot_alt)\n",
    "Plots.savefig(erg_plot_alt, string(\".\\\\\", model_instance, \"\\\\ergodicity_convergence_all_alt.eps\"));\n",
    "\n",
    "display(r_plot)\n",
    "Plots.savefig(r_plot, string(\".\\\\\", model_instance, \"\\\\ergodicity_convergence.eps\"));\n",
    "Plots.savefig(rho_plot, string(\".\\\\\", model_instance, \"\\\\ergodicity_convergence_rho.eps\"));\n",
    "Plots.savefig(G_plot, string(\".\\\\\", model_instance, \"\\\\ergodicity_convergence_G.eps\"));\n",
    "\n",
    "display(rho_plot)\n",
    "display(G_plot)\n",
    "display(P_plot)\n",
    "#Plots.plot(rho_plot, G_plot, P_plot)\n",
    "#Plots.savefig(string(\".\\\\\", model_instance, \"\\\\ergodicity_convergence_decomposed.pdf\"));\n",
    "\n",
    "display(MU_plot)\n",
    "Plots.plot(MU_plot)\n",
    "Plots.savefig(string(\".\\\\\", model_instance, \"\\\\ergodicity_payoff.eps\"));\n",
    "\n",
    "display(MD_plot)\n",
    "Plots.plot(MD_plot)\n",
    "Plots.savefig(string(\".\\\\\", model_instance, \"\\\\ergodicity_damage.eps\"));\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# End of task: plot the Monte Carlo runs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Task: initialise some required packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Pkg;\n",
    "Pkg.build(\"Plots\")\n",
    "Pkg.build(\"GR\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import Pkg;\n",
    "\n",
    "Pkg.add(\"Plots\"); Pkg.add(\"PyPlot\"); Pkg.add(\"LaTeXStrings\"); Pkg.add(\"Colors\"); Pkg.add(\"GR\");\n",
    "Pkg.add(\"ColorSchemes\")\n",
    "Pkg.add(\"JuMP\"); Pkg.add(\"Ipopt\");\n",
    "Pkg.add(\"BasisMatrices\"); \n",
    "Pkg.add(\"Parameters\");\n",
    "Pkg.add(\"JLD2\"); Pkg.add(\"FileIO\"); Pkg.add(\"CSV\"); Pkg.add(\"DataFrames\")\n",
    "Pkg.add(\"StatsBase\")\n",
    "Pkg.add(\"GRUtils\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# End of task: initialise required packages."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.0.5",
   "language": "julia",
   "name": "julia-1.0"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.0.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
